# R Script for Teargas and Selfie Cams
# Last updated: December 29, 2020


## EXPLANATION: LINES 7 TO 164 LOAD THE DATA FOR STUDY 1 IN THE PAPER. ##

##################################################################
##################################################################
###################### STUDY 1 - 2016 ############################
##################################################################
##################################################################

#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
rm(list=ls())
#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ipak <- function(pkg){  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}

packages <- c("memisc", "car", "stm", "tm", "wordcloud", "stargazer", "dplyr", "stminsights")
ipak(packages)

library(logr)
library(car)
library(stm)
library(stminsights)
library(stargazer)
library(dplyr)
library(ggplot2)

# log for output
TandSC_log <- file.path(tempdir(), "GKF.log")
logfile <- log_open(TandSC_log)
#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# setwd() -- SET TO APPROPRIATE FOLDER

#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
##LOADING THE EXPERIMENTAL DATA

exp_2016 <- read.csv("FGK_51916.csv")
head(exp_2016)
dim(exp_2016)
#removing the first row which has question text
exp_2016 <- exp_2016[-1,]
exp_2016 <- exp_2016[exp_2016$consent==1,]
exp_2016[exp_2016$consent==2,]

##CREATING VARIABLES

##Creating one column for all open ended responses
exp_2016$all_text <- paste(exp_2016$opentextfavor, exp_2016$opentextneither, exp_2016$opentextoppose)

##coding the violence variable, pooling together all violent conditions
exp_2016$violence <- ifelse(exp_2016$violentEgypt==1 | exp_2016$violentUkraine==1 | exp_2016$violentHK == 1 | exp_2016$Q143 == 1 | exp_2016$Q145 == 1 | exp_2016$Q147 == 1 | exp_2016$Q150 == 1 | exp_2016$Q152 == 1 | exp_2016$Q154 == 1 | exp_2016$Q157==1 | exp_2016$Q159==1 | exp_2016$Q161 ==1 | exp_2016$Q164==1 | exp_2016$Q166 == 1 | exp_2016$Q168 ==1 | exp_2016$Q171==1 | exp_2016$Q173 == 1 | exp_2016$Q175 == 1 | exp_2016$Q178==1 | exp_2016$Q180 ==1 | exp_2016$Q182 == 1, 1, 0)
sum(exp_2016$violence)

###code the policy variable

exp_2016$policy2 <- as.numeric(as.character(exp_2016$policy))
exp_2016$favor_sanctions <- recode_factor(exp_2016$policy2, '1'="7", '3'="6", '4'="5", '5'="4", '6'="3", '7'="2", '8'="1", .missing="NA")
#exp_2016$favor_sanctions <- recode(exp_2016$policy2, "1='7'; 3='6'; 4='5'; 5='4'; 6='3'; 7='2'; 8='1', 'NA'=NA", as.factor.result=FALSE)

##binary for comparability
exp_2016$favor_sanctions <- as.numeric(as.character(exp_2016$favor_sanctions))
exp_2016$support_sanctions <- ifelse(exp_2016$favor_sanctions>4, 1, 0)

##converting levels into numeric values
exp_2016$education <- as.numeric(as.character(exp_2016$edu))
exp_2016$college_degree <- ifelse(exp_2016$education>=5, 1, 0)

exp_2016$income.num <- as.numeric(as.character(exp_2016$income))
exp_2016$mid_income <- ifelse(exp_2016$income.num >= 5, 1, 0)

exp_2016$news.num <- as.numeric(as.character(exp_2016$news))

#make an age variable
exp_2016$age <- 2016-as.numeric(as.character(exp_2016$birthY))
#gender variable
exp_2016$male <- ifelse(exp_2016$gender==1, 1, 0)

#make ethnicity variable
exp_2016$white <- ifelse(exp_2016$ethnicity==2, 1, 0)

##make a voting variable
exp_2016$vote.num <- ifelse(exp_2016$vote==2, 0, 1)

##make partisanship variables
exp_2016$democrat <- ifelse(exp_2016$partyID == 1 | exp_2016$leaning == 2, 1, 0)
exp_2016$republican <-ifelse(exp_2016$democrat==1, 0, 1)

##make a protester variable
exp_2016$protester <- ifelse(exp_2016$protest_1 == 1 | exp_2016$protest_2 == 1 | exp_2016$protest_3 == 1 | exp_2016$protest_4 == 1 | exp_2016$protest_5 == 1 | exp_2016$protest_6 == 1 | exp_2016$protest_7 == 1, 1, 0)

#affinity towards regions
exp_2016$Asia_thermometer <- as.numeric(as.character(exp_2016$thermometer_4))
exp_2016$ME_thermometer <- as.numeric(as.character(exp_2016$thermometer_6))
exp_2016$EE_thermometer<- as.numeric(as.character(exp_2016$thermometer_2))

#code which protest they saw
exp_2016$egypt <- ifelse(exp_2016$nonviolentEgypt == 1 | exp_2016$Q142 == 1 | exp_2016$Q149 == 1 | exp_2016$Q156 == 1 | exp_2016$Q163 == 1 | exp_2016$Q170 == 1 | exp_2016$Q177 == 1 | exp_2016$violentEgypt == 1 | exp_2016$Q143 == 1 | exp_2016$Q150 == 1 | exp_2016$Q157 == 1 | exp_2016$Q164 == 1 | exp_2016$Q171 == 1 | exp_2016$Q178 == 1, 1, 0)

exp_2016$ukraine <- ifelse(exp_2016$nonviolentUkraine == 1 | exp_2016$Q144 == 1 | exp_2016$Q151 == 1 | exp_2016$Q158 == 1 | exp_2016$Q165 == 1 | exp_2016$Q172 == 1 | exp_2016$Q179 == 1 | exp_2016$violentUkraine == 1 | exp_2016$Q145 == 1 | exp_2016$Q152 == 1 | exp_2016$Q159 == 1 | exp_2016$Q166 == 1 | exp_2016$Q173 == 1 | exp_2016$Q180 == 1, 1, 0)

exp_2016$hk <- ifelse(exp_2016$nonviolentHK == 1 | exp_2016$Q146 == 1 | exp_2016$Q153 == 1 | exp_2016$Q160 == 1 | exp_2016$Q167 == 1 | exp_2016$Q174 == 1 | exp_2016$Q181 == 1 | exp_2016$violentHK == 1 | exp_2016$Q147 == 1 | exp_2016$Q154 == 1 | exp_2016$Q161 == 1 | exp_2016$Q168 == 1 | exp_2016$Q175 == 1 | exp_2016$Q182 == 1, 1, 0)

#create one variable for feeling thermometer towards region viewed
exp_2016$thermometer <- NULL

for(i in 1:nrow(exp_2016)){
  if(exp_2016$ukraine[i]==1) {
    exp_2016$thermometer[i] <- exp_2016$EE_thermometer[i]
  } else if (exp_2016$egypt[i]==1) {
    exp_2016$thermometer[i] <- exp_2016$ME_thermometer[i]
  } else {
    exp_2016$thermometer[i] <- exp_2016$Asia_thermometer[i]
  }
}

##make a money donation variable
exp_2016$donation.1.num <- as.numeric(as.character(exp_2016$donation.1))
exp_2016$donation.2.num <- as.numeric(as.character(exp_2016$donation.2))

##recoding NAs as 0 so that I can collapse the two vectors into one
exp_2016$donation.1.num[is.na(exp_2016$donation.1.num)] <- 0
exp_2016$donation.2.num[is.na(exp_2016$donation.2.num)] <- 0
exp_2016$money <- exp_2016$donation.1.num + exp_2016$donation.2.num
table(exp_2016$money)


##make a variable for behavior, willing to take action on behalf of protesters, agnostic on costliness of action 
exp_2016$leader <- ifelse(exp_2016$action.1_1 == 1 | exp_2016$action.2_7 == 1, 1, 0)
exp_2016$organize <- ifelse(exp_2016$action.1_2==1 | exp_2016$action.2_6==1, 1, 0)
exp_2016$protest <- ifelse(exp_2016$action.1_3==1 | exp_2016$action.2_4==1, 1,0 )
exp_2016$petition <- ifelse(exp_2016$action.1_7==1 | exp_2016$action.2_1==1,1,0 )
exp_2016$letter <- ifelse(exp_2016$action.1_4==1 | exp_2016$action.2_3==1,1,0)
exp_2016$dmoney <- ifelse(exp_2016$action.1_6==1 | exp_2016$action.2_2==1, 1,0)
exp_2016$none <- ifelse(exp_2016$action.1_5==1 | exp_2016$action.2_5 == 1,1,0)

#make an index of political behavior
exp_2016$index <- NULL
for(i in 1:nrow(exp_2016))
{
  if (exp_2016$leader[i]==1) {
    exp_2016$index[i] <- 6
  } else if (exp_2016$organize[i]==1) {
    exp_2016$index[i] <- 5
  } else if (exp_2016$protest[i]==1) {
    exp_2016$index[i] <- 4
  } else if (exp_2016$petition[i]==1) {
    exp_2016$index[i] <- 3
  } else if (exp_2016$letter[i]==1) {
    exp_2016$index[i] <- 2
  } else if (exp_2016$dmoney[i]==1) {
    exp_2016$index[i] <- 1
  } else {
    exp_2016$index[i] <- 0
  }
}
table(exp_2016$index)

# create a variable for if respondent would take ANY action on behalf of protesters, based on "none" category
exp_2016$pol_behavior <- ifelse(exp_2016$action.1_5==1 | exp_2016$action.2_5 == 1,0,1)
table(exp_2016$pol_behavior)

#########################Analysis##################################

## EXPLANATION: LINES 178 TO 474 GIVE RESULTS FOR STUDY 1 IN THE PAPER, STM STARTS AT 368 ##
## RUN STUDY 1 CODE IN FULL (LINES 7 TO 870) TO MAKE SURE EVERYTHING WORKS AS CODE OFTEN RELIES ON CODE BEFORE IT ##

## EXPLANATION: THE MEANS IN THE PAPER (WRITTEN IN TEXT ON PAGE 6) CAN BE FOUND IN LINES 187 to 193 ##

## EXPLANATION: LINES 262 TO 288 GIVES PAPER FIGURE 1 PLOT ### 

## EXPLANATION: LINES 255 TO 257 GIVE APPENDIX TABLE 4 (STUDY 1 REGRESSION TABLE) ##


## T-TESTS AND REGRESSIONS

## Support for Sanctions
  log_print("Study 1, support for sanctions, difference in means")
    
##check difference in means for effect of violence on sanctions support
diff.in.means_95 <- t.test(exp_2016[exp_2016$violence==1,]$support_sanctions, exp_2016[exp_2016$violence==0,]$support_sanctions)
diff.in.means_90 <- t.test(exp_2016[exp_2016$violence==1,]$support_sanctions, exp_2016[exp_2016$violence==0,]$support_sanctions, conf.level = 0.9)

# mean for nonviolent videos

diff.in.means_95$estimate[2]

# mean for violent videos

diff.in.means_95$estimate[1]

effect1 <- diff.in.means_95$estimate[1]-diff.in.means_95$estimate[2]

log_print(diff.in.means_95)

##bivariate regression with violence as IV 
H1reg.binary.sanc <- lm(support_sanctions~violence, data=exp_2016)
summary(H1reg.binary.sanc)

##regression controlling for region and violence
H1reg.region.binary.sanc <- lm(support_sanctions~violence + as.factor(ukraine) + as.factor(egypt), data=exp_2016)
summary(H1reg.region.binary.sanc)

##regression with all covariates
H1reg.full.binary.sanc <- lm(support_sanctions~violence + as.factor(ukraine) + as.factor(egypt) +  age + male + white + vote.num + republican + college_degree + mid_income  + protester + news.num  + thermometer, data=exp_2016)
summary(H1reg.full.binary.sanc)

## Donations
##check difference in means for effect of violence on donations
diff.in.means2_95 <- t.test(exp_2016[exp_2016$violence==1,]$money, exp_2016[exp_2016$violence==0,]$money)
diff.in.means2_90 <- t.test(exp_2016[exp_2016$violence==1,]$money, exp_2016[exp_2016$violence==0,]$money, conf.level = 0.9)

effect2 <- diff.in.means2_95$estimate[1]-diff.in.means2_95$estimate[2]

##bivariate regression with violence as IV 
H1reg.binary.don <- lm(money~violence, data=exp_2016)
summary(H1reg.binary.don)

##regression controlling for region and violence
H1reg.region.binary.don <- lm(money~violence + as.factor(ukraine) + as.factor(egypt), data=exp_2016)
summary(H1reg.region.binary.don)

##regression with all covariates
H1reg.full.binary.don <- lm(money~violence + as.factor(ukraine) + as.factor(egypt) +  age + male + white + vote.num + republican + college_degree + mid_income  + protester + news.num  + thermometer, data=exp_2016)
summary(H1reg.full.binary.don)

## Political Behavior
##check difference in means for effect of violence on political behavior
diff.in.means3_95 <- t.test(exp_2016[exp_2016$violence==1,]$pol_behavior, exp_2016[exp_2016$violence==0,]$pol_behavior)
diff.in.means3_90 <- t.test(exp_2016[exp_2016$violence==1,]$pol_behavior, exp_2016[exp_2016$violence==0,]$pol_behavior, conf.level = 0.9)

effect3 <- diff.in.means3_95$estimate[1]-diff.in.means3_95$estimate[2]

##bivariate regression with violence as IV 
H1reg.binary.behav <- lm(pol_behavior~violence, data=exp_2016)
summary(H1reg.binary.behav)

##regression controlling for region and violence
H1reg.region.binary.behav <- lm(pol_behavior~violence + as.factor(ukraine) + as.factor(egypt), data=exp_2016)
summary(H1reg.region.binary.behav)

##regression with all covariates
H1reg.full.binary.behav <- lm(pol_behavior~violence + as.factor(ukraine) + as.factor(egypt) +  age + male + white + vote.num + republican + college_degree + mid_income  + protester + news.num  + thermometer, data=exp_2016)
summary(H1reg.full.binary.behav)


#bivariate regression
summary(H1reg.binary.sanc)
summary(H1reg.binary.don)
summary(H1reg.binary.behav)

#Generate table for appendix
## CODE FOR APPENDIX TABLE 4
stargazer(H1reg.binary.sanc,  H1reg.region.binary.sanc, H1reg.full.binary.sanc, H1reg.binary.don,  H1reg.region.binary.don, H1reg.full.binary.don, H1reg.binary.behav,  H1reg.region.binary.behav, H1reg.full.binary.behav, digits=3, single.row=TRUE, title="Study 1: Effects of Violent Police Response")

log_print("Study 1 Regression Table, code for Appendix Table 4")
log_print(stargazer(H1reg.binary.sanc,  H1reg.region.binary.sanc, H1reg.full.binary.sanc, H1reg.binary.don,  H1reg.region.binary.don, H1reg.full.binary.don, H1reg.binary.behav,  H1reg.region.binary.behav, H1reg.full.binary.behav, digits=3, single.row=TRUE, title="Study 1: Effects of Violent Police Response"))

## PLOT FOR FIGURE 1
###plot violence effects for all three DVs
dev.off()
par(mfrow=c(3,1))
plot(effect1,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_95$conf.int[1],3,diff.in.means_95$conf.int[2],3,lwd=2)
segments(diff.in.means_90$conf.int[1],3,diff.in.means_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Support for Sanctions",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect2,3,cex=3,pch=20,xlim=c(-15,15),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(diff.in.means2_95$conf.int[1],3,diff.in.means2_95$conf.int[2],3,lwd=2)
segments(diff.in.means2_90$conf.int[1],3,diff.in.means2_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Money Donations",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-15,15,by=5)),labels=c("-15", "-10", "-5", "0", "5", "10", "15"))
mtext("Dollar Amount Change",side=1,cex=1.1,line=2.3)

plot(effect3,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(diff.in.means3_95$conf.int[1],3,diff.in.means3_95$conf.int[2],3,lwd=2)
segments(diff.in.means3_90$conf.int[1],3,diff.in.means3_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Political Behavior",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

###check effects when DVs are ordinal
#not applicable to money donation (cont.) 

# support for sanctions

##regression controlling for region and violence
H1reg.ordinal.sanc <- lm(favor_sanctions~violence, data=exp_2016)
summary(H1reg.ordinal.sanc)

##regression controlling for region and violence
H1reg.region.ordinal.sanc <- lm(favor_sanctions~violence + as.factor(ukraine) + as.factor(egypt), data=exp_2016)
summary(H1reg.region.ordinal.sanc)

##regression with all covariates
H1reg.full.ordinal.sanc <- lm(favor_sanctions~violence + as.factor(ukraine) + as.factor(egypt) +  age + male + white + vote.num + republican + college_degree + mid_income  + protester + news.num  + thermometer, data=exp_2016)
summary(H1reg.full.ordinal.sanc)

# political behavior

##regression controlling for region and violence
H1reg.ordinal.behav <- lm(index~violence, data=exp_2016)
summary(H1reg.ordinal.behav)

##regression controlling for region and violence
H1reg.region.ordinal.behav <- lm(index~violence + as.factor(ukraine) + as.factor(egypt), data=exp_2016)
summary(H1reg.region.ordinal.behav)

##regression with all covariates
H1reg.full.ordinal.behav <- lm(index~violence + as.factor(ukraine) + as.factor(egypt) +  age + male + white + vote.num + republican + college_degree + mid_income  + protester + news.num  + thermometer, data=exp_2016)
summary(H1reg.full.ordinal.behav)

###we can also change effects to results of binary IV regressions with covariates in plots
effect1a <- summary(H1reg.full.binary.sanc)$coef[2]
effect2a <- summary(H1reg.full.binary.don)$coef[2]
effect3a <- summary(H1reg.full.binary.behav)$coef[2]

ui1_95 <- effect1a + 1.96 * summary(H1reg.full.binary.sanc)$coefficients[,2][2]
li1_95 <- effect1a - 1.96 * summary(H1reg.full.binary.sanc)$coefficients[,2][2]
ui2_95 <- effect2a + 1.96 * summary(H1reg.full.binary.don)$coefficients[,2][2]
li2_95 <- effect2a - 1.96 * summary(H1reg.full.binary.don)$coefficients[,2][2]
ui3_95 <- effect3a + 1.96 * summary(H1reg.full.binary.behav)$coefficients[,2][2]
li3_95 <- effect3a - 1.96 * summary(H1reg.full.binary.behav)$coefficients[,2][2]

ui1_90 <- effect1a + 1.645 * summary(H1reg.full.binary.sanc)$coefficients[,2][2]
li1_90 <- effect1a - 1.645 * summary(H1reg.full.binary.sanc)$coefficients[,2][2]
ui2_90 <- effect2a + 1.645 * summary(H1reg.full.binary.don)$coefficients[,2][2]
li2_90 <- effect2a - 1.645 * summary(H1reg.full.binary.don)$coefficients[,2][2]
ui3_90 <- effect3a + 1.645 * summary(H1reg.full.binary.behav)$coefficients[,2][2]
li3_90 <- effect3a - 1.645 * summary(H1reg.full.binary.behav)$coefficients[,2][2]


###plot violence effects using regressions with covariates

par(mfrow=c(3,1))
plot(effect1a,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(li1_95,3,ui1_95,3,lwd=2)
segments(li1_90,3,ui1_90,3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Support for Sanctions",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect2a,3,cex=3,pch=20,xlim=c(-15,15),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li2_95,3,ui2_95,3,lwd=2)
segments(li2_90,3,ui2_90,3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Money Donations",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-15,15,by=5)),labels=c("-15", "-10", "-5", "0", "5", "10", "15"))
mtext("Dollar Amount Change",side=1,cex=1.1,line=2.3)

plot(effect3a,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li3_95,3,ui3_95,3,lwd=2)
segments(li3_90,3,ui3_90,3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Political Behavior",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

#########################Structural Topic Model##################################

## STRUCTURAL TOPIC MODEL

## EXPLANATION: LINES 444 TO 461 GIVE PAPER FIGURE 2 (STM TOPICS) ##

## EXPLANATION: LINES 466 TO 474 GIVE PAPER FIGURE 3 (STM EFFECTS PLOT)


set.seed(02138)

##processing the data - lower case, removes punctuation, stemming etc.

temp<-textProcessor(documents=exp_2016$all_text,metadata=exp_2016, stem=TRUE, removestopwords=TRUE, language="en")

##preparing an object for stm analysis
out <- prepDocuments(temp$documents, temp$vocab, temp$meta)

documents <- out$documents

vocab <- out$vocab

meta <- out$meta


## run a for loop to pick number of topics (k) for the model with treatment status as covariate
## produces plots with frex words with topic number 3 - 10
##stm documentation recommends 3-10 as a useful starting range for survey experiments, 3 minimum, more than 10 too hard to interpret


for (i in 3:10) {
  temp <- stm(documents, vocab, K=i, prevalence=~violence, data=meta, seed=02138)
  pdf(file = paste("Frex.Words_replicate.",i,".pdf",collapse="",sep="")) 
  par(mar=c(0,0,0,0))
  plot.STM(temp,type="labels", n=10, labeltype="frex")
  dev.off()
}

##choose three topics after examining output because after > 3, topics hard to distinguish

##model selection for k=3

mod.select <- selectModel(out$documents, out$vocab, K=3, prevalence=~violence, data=meta, max.em.its=500, runs=50, seed=02138)

par(mfrow=c(1,1))
plotModels(mod.select) ##choose model on the upper right (based exclusivity and coherence)

violenceModel <- mod.select$runout[[4]] ##choose the fourth model

log_print("Study 1, structural topic model")
log_print(violenceModel)

##examine output of the model, words associated with each topic to infer topic labels

labelTopics(violenceModel, n=20)

##check representative documents for each topic as a validation step for labelling topics

##topic 1
thoughts1 <- findThoughts(violenceModel,texts=meta$all_text, topics=1, n=10)$docs[[1]]
plotQuote(thoughts1[1], text.cex=2)

##topic 2
thoughts2 <- findThoughts(violenceModel,texts=meta$all_text, topics=2, n=10)$docs[[1]]
plotQuote(thoughts2[3], text.cex=2)

##topic 3
thoughts3 <- findThoughts(violenceModel,texts=meta$all_text, topics=3, n=10)$docs[[1]]
plotQuote(thoughts3[2], text.cex=2)

log_print("findThoughts for each of the STM three topics")

log_print(thoughts1)
log_print(thoughts2)
log_print(thoughts3)

## PLOT FOR FIGURE 2
###plot frequent and exclusive words for each topic 
plot.STM(violenceModel,type="labels", n=20, labeltype="frex", topic.names=c("Topic 1: Non-Interventionism", "Topic 2: Don't Know/Not enough info", "Topic 3: Human Rights"))


############# estimate effect of police violence on topic prevalence ############

violModel_effect <- estimateEffect(1:3 ~ violence, violenceModel, meta)
summary(violModel_effect)
violModel_effect$varlist[1]
get_effects(violModel_effect, variable= violModel_effect$varlist[1], cov_val1 = 0, cov_val2 = 1, type= "difference")

stm_90 <- get_effects(violModel_effect, variable= violModel_effect$varlist[1], cov_val1 = 0, cov_val2 = 1, type= "difference", ci = 0.9)

log_print("STM model outputs")

log_print(summary(violModel_effect))
log_print(stm_90)

## plot the effect of bombing on proportion of documents dedicated in each topic; confidence intervals 95% 


## PLOT FOR FIGURE 3

par(mfrow=c(1,1))
plot.estimateEffect(violModel_effect, covariate = "violence", topics = 1:3, model=violenceModel, method="difference", ci.level=0.95,
                    cov.value1="1",cov.value2="0", main="", xlab="Difference in Topic Proportions between Violent and Nonviolent Police Response",
                    xlim=c(-.2,.2), labeltype = "custom", custom.labels = c("Topic 1: Non-Interventionism", "Topic 2: Don't Know/Not enough info", "Topic 3: Human Rights"))
segments(as.numeric(-stm_90[1,3]),3,as.numeric(-stm_90[1,4]),3,lwd=4) # adding 90% CI, Topic 1
segments(as.numeric(-stm_90[2,3]),2,as.numeric(-stm_90[2,4]),2,lwd=4) # adding 90% CI, Topic 2
segments(as.numeric(-stm_90[3,3]),1,as.numeric(-stm_90[3,4]),1,lwd=4) # adding 90% CI, Topic 3


#########################Checks##################################

## MANIPULATION CHECK, ATTENTION CHECK, BALANCE CHECKS AND ROBUSTNESS CHECKS

## EXPLANATION: LINES 503 TO 568 GIVE APPENDIX FIGURE 1 TOP RIGHT ##

## EXPLANATION: LINES 569 TO 598 GIVE APPENDIX FIGURE 1 TOP LEFT ##

## EXPLANATION: LINES 602 TO 638 GIVE APPENDIX FIGURE 1 BOTTOM ##

## EXPLANATION: LINES 707 TO 753 GIVE DATA FOR APPENDIX TABLE 2 ##

## EXPLANATION: LINES 667 TO 705 GIVE DATA FOR APPENDIX TABLE 3 ##

## EXPLANATION: LINES 837 TO 888 GIVE APPENDIX FIGURE 2 ##

## EXPLANATION: LINES 757 TO 775 GIVE APPENDIX ROBUSTNESS CHECK ESTIMATE OF ~27% (WRITTEN IN TEXT ON APPENDIX PAGE 3) ##

## EXPLANATION: LINES 656 TO 661 GIVE APPENDIX ATTENTION CHECK ESTIMATE OF ~97% (WRITTEN IN TEXT ON APPENDIX PAGE 2) ##


#manipulation check using pre-test data
#respondents were asked about levels of violence in each of the 6 videos

# setwd() -- SET TO APPROPRIATE FOLDER

#load data
qualtrics <- read.csv("PreTest42116.csv")


names(data)

qualtrics$me <- ifelse(qualtrics$MEpeaceful==1 | qualtrics$MEviolent==1 | qualtrics$NonViolentME==1 | qualtrics$violentME==1, 1, 0)
qualtrics$me[is.na(qualtrics$me)] <- 0

qualtrics$vio <- ifelse(qualtrics$violentME==1 | qualtrics$violentUkraine==1 | qualtrics$violentHK==1, 1, 0)
qualtrics$vio[is.na(qualtrics$vio)] <- 0



qualtrics$nv <- ifelse(qualtrics$NonViolentME==1 | qualtrics$NonViolentUkraine==1 | qualtrics$NonViolentHK==1, 1, 0)
qualtrics$nv[is.na(qualtrics$nv)] <- 0


##################police violence
library(ggplot2)

#creating plots to determine whether there is a difference in perceived levels of police violence in videos

##police violence - violent group
polvioVMEa <- na.omit(as.numeric(as.character(qualtrics$polvioVME_1)))
polvioVMEa
mean(polvioVMEa)
polvioVHKa <- na.omit(as.numeric(as.character(qualtrics$polvioVHK_1)))
mean(polvioVHKa)
polvioVukraineA <- na.omit(as.numeric(as.character(qualtrics$polvioVukraine_1)))
mean(polvioVukraineA)

vhkvuktt <- t.test(polvioVHKa, polvioVukraineA, mu=0) ##do not reject the null for any of them
#point
vhkvuktt.point <- vhkvuktt$estimate[1]-vhkvuktt$estimate[2]

#lower bound
vhkvuktt$conf.int[1]
#uper bound
vhkvuktt$conf.int[2]

vhkvmett <- t.test(polvioVHKa, polvioVMEa)
vhkvmett.point <- vhkvmett$estimate[1] - vhkvmett$estimate[2]

vukvmett <- t.test(polvioVukraineA, polvioVMEa)
vukvmett.point <- vukvmett$estimate[1] - vukvmett$estimate[2]
as.numeric(vukvmett.point)

points <- c(as.numeric(vhkvuktt.point), as.numeric(vhkvmett.point), as.numeric(vukvmett.point))
lower95 <- c(vhkvuktt$conf.int[1], vhkvmett$conf.int[1], vukvmett$conf.int[1])
upper95 <- c(vhkvuktt$conf.int[2], vhkvmett$conf.int[2], vukvmett$conf.int[2])
x <- c(1:3)
df <- as.data.frame(cbind(points, lower95, upper95, x))
class(df)

df
p = qplot(factor(x=1:nrow(df)), y=points, data=df) 

p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95))
p = p + ylab("Difference-in-means") + ggtitle("Levels of Police Violence in Videos: Violent Group")
p = p + scale_x_discrete(breaks = 1:3, labels=c("Asia-Eur","Asia-ME","Eur-ME")) + xlab(NULL)
p = p + geom_hline(yintercept=0, linetype="dotted")
p

## PLOT FOR APPENDIX FIGURE 1 - first plot
print(p)


##nonviolent group
polvioNVMEa <- na.omit(as.numeric(as.character(qualtrics$polvioNVME_1)))
polvioNVHKa <- na.omit(as.numeric(as.character(qualtrics$polvioNVHK_1)))
polvioNVukraineA <- na.omit(as.numeric(as.character(qualtrics$polvioNVukraine_1)))

nvhkuk <- t.test(polvioNVHKa, polvioNVukraineA)
nvhkuk.point <- nvhkuk$estimate[1] - nvhkuk$estimate[2]
nvhkmett <- t.test(polvioNVHKa, polvioNVMEa)
nvhkmett.point <- nvhkmett$estimate[1] - nvhkmett$estimate[2]
nvuknvmett <- t.test(polvioNVukraineA, polvioNVMEa)
nvuknvmett.point <- nvuknvmett$estimate[1] - nvuknvmett$estimate[2]

points2 <- c(as.numeric(nvhkuk.point), as.numeric(nvhkmett.point), as.numeric(nvuknvmett.point))
lower952 <- c(nvhkuk$conf.int[1], nvhkmett$conf.int[1], nvuknvmett$conf.int[1])
upper952 <- c(nvhkuk$conf.int[2], nvhkmett$conf.int[2], nvuknvmett$conf.int[2])
x <- c(1:3)
df2 <- as.data.frame(cbind(points2, lower952, upper952, x))

df2
p2 = qplot(factor(x=1:nrow(df2)), y=points2, data=df2) 

p2 = p2 + geom_pointrange(aes(ymin=lower952,ymax=upper952))
p2 = p2 + ylab("Difference-in-means") + ggtitle("Levels of Police Violence in Videos: Non-Violent Group")
p2 = p2 + scale_x_discrete(breaks = 1:3, labels=c("Asia-Eur","Asia-ME","Eur-ME")) + xlab(NULL)
p2 = p2 + geom_hline(yintercept=0, linetype="dotted")

## PLOT FOR APPENDIX FIGURE 1 - second plot
p2



##ukraine
polvioUkraineNVa <- na.omit(as.numeric(as.character(qualtrics$polvioUkraineNV_1)))
polvioUkraineVa <- na.omit(as.numeric(as.character(qualtrics$polvioUkraineV_1)))

ukraine <- t.test(polvioUkraineVa, polvioUkraineNVa) ##reject the null
ukraine.point <- as.numeric(ukraine$estimate[1] - ukraine$estimate[2]) 

##middle east
polvioMEVa <- na.omit(as.numeric(as.character(qualtrics$polvioMEV_1))) 
polvioMENVa <- na.omit(as.numeric(as.character(qualtrics$polvioMENV_1))) 

me <- t.test(polvioMEVa, polvioMENVa) ###reject the null
me.point <- as.numeric(me$estimate[1] - me$estimate[2])

##hong kong
polvioHKNVa <- na.omit(as.numeric(as.character(qualtrics$polvioHKNV_1)))
polvioHKVa <- na.omit(as.numeric(as.character(qualtrics$polvioHKviolent_1)))

hk <- t.test(polvioHKVa, polvioHKNVa) ##reject the null
hk.point <- as.numeric(hk$estimate[1] - hk$estimate[2])

points3 <- c(ukraine.point, me.point, hk.point)
lower953 <- c(ukraine$conf.int[1], me$conf.int[1], hk$conf.int[1])
upper953 <- c(ukraine$conf.int[2], me$conf.int[2], hk$conf.int[2])
x <- c(1:3)
df3 <- as.data.frame(cbind(points3, lower953, upper953, x))

df3
p3 = qplot(factor(x=1:nrow(df3)), y=points3, data=df3, ylim=c(-1,6)) 

p3 = p3 + geom_pointrange(aes(ymin=lower953,ymax=upper953))
p3 = p3 + ylab("Difference-in-means") + ggtitle("Levels of Police Violence in Videos: Within Country")
p3 = p3 + scale_x_discrete(breaks = 1:3, labels=c("Violent Eur-NonViolent Eur","Violent ME-NonViolent ME","Violent Asia-NonViolent Asia")) + xlab(NULL)
p3 = p3 + geom_hline(yintercept=0, linetype="dotted")

## PLOT FOR APPENDIX FIGURE 1 - plot 3
p3


###each video - police violence levels

##middle east nonviolent
mean(c(polvioNVMEa, polvioMENVa))
###middle east violent
mean(c(polvioMEVa, polvioMEVa))
##ukraine nonviolent
mean(c(polvioUkraineNVa, polvioNVukraineA))
##ukraine violent
mean(c(polvioUkraineVa, polvioVukraineA))
##hong kong nonviolent
mean(c(polvioHKNVa, polvioNVHKa))
##hong kong violent
mean(c(polvioHKVa, polvioVHKa))

# attention check

table(exp_2016$attentioncheck)/nrow(exp_2016)
# 97% of respondents got the question right
log_print("Study 1, attention check")
log_print(table(exp_2016$attentioncheck)/nrow(exp_2016))

############Balance Tests#################

# balance checks on control variables

HK <- exp_2016[exp_2016$hk == 1,]
Egypt <- exp_2016[exp_2016$egypt == 1,]
Ukraine <- exp_2016[exp_2016$ukraine == 1,]

t.test(HK$age)
t.test(Egypt$age)
t.test(Ukraine$age)

t.test(HK$male)
t.test(Egypt$male)
t.test(Ukraine$male)

t.test(HK$white)
t.test(Egypt$white)
t.test(Ukraine$white)

t.test(HK$vote.num)
t.test(Egypt$vote.num)
t.test(Ukraine$vote.num)

t.test(HK$republican)
t.test(Egypt$republican)
t.test(Ukraine$republican)

t.test(HK$college_degree)
t.test(Egypt$college_degree)
t.test(Ukraine$college_degree)

t.test(HK$mid_income)
t.test(Egypt$mid_income)
t.test(Ukraine$mid_income)

t.test(HK$protester)
t.test(Egypt$protester)
t.test(Ukraine$protester)

t.test(HK$news.num)
t.test(Egypt$news.num)
t.test(Ukraine$news.num)

violent <- exp_2016[exp_2016$violence == 1,]
nonviolent <- exp_2016[exp_2016$violence == 0,]

t.test(violent$age, nonviolent$age)
sd(violent$age)
sd(nonviolent$age)
# no significant difference

t.test(violent$male, nonviolent$male)
sd(violent$male)
sd(nonviolent$male)
# no significant difference

t.test(violent$white, nonviolent$white)
sd(violent$white)
sd(nonviolent$white)
# no significant difference

t.test(violent$vote.num, nonviolent$vote.num)
sd(violent$vote.num)
sd(nonviolent$vote.num)
# no significant difference

t.test(violent$republican, nonviolent$republican)
sd(violent$republican)
sd(nonviolent$republican)
# no significant difference

t.test(violent$college_degree, nonviolent$college_degree)
sd(violent$college_degree)
sd(nonviolent$college_degree)
# no significant difference

t.test(violent$mid_income, nonviolent$mid_income)
sd(violent$mid_income)
sd(nonviolent$mid_income)
# no significant difference

t.test(violent$protester, nonviolent$protester)
sd(violent$protester)
sd(nonviolent$protester)
# no significant difference

t.test(violent$news.num, nonviolent$news.num)
sd(violent$news.num)
sd(nonviolent$news.num)
# no significant difference

############Robustness Checks#################

##Robustness check - analysis excluding respondents who correctly identified the location of the protest 
#checking on what respondents wrote when asked about the protest region
head(exp_2016$locationcheck)
unique(exp_2016$locationcheck)
table(grepl("egypt", exp_2016$locationcheck, ignore.case = TRUE))
table(grepl("china", exp_2016$locationcheck, ignore.case = TRUE))
table(grepl("hong kong", exp_2016$locationcheck, ignore.case = TRUE))
table(grepl("hk", exp_2016$locationcheck, ignore.case = TRUE))
table(grepl("ukraine", exp_2016$locationcheck, ignore.case = TRUE))
table(grepl("cairo", exp_2016$locationcheck, ignore.case = TRUE))
table(grepl("kiev", exp_2016$locationcheck, ignore.case = TRUE))
table(grepl("eqypt", exp_2016$locationcheck, ignore.case = TRUE))
## so 50 egypt,1 eqypt, 133china, 38 hong kong, 55 ukraine, 6 cairo, 4kiev

exp_2016$awarelocation <- ifelse(grepl("egypt", exp_2016$locationcheck, ignore.case = TRUE) | grepl("china", exp_2016$locationcheck, ignore.case = TRUE) | grepl("hong kong", exp_2016$locationcheck, ignore.case = TRUE) | grepl("ukraine", exp_2016$locationcheck, ignore.case = TRUE) | grepl("cairo", exp_2016$locationcheck, ignore.case = TRUE) | grepl("kiev", exp_2016$locationcheck, ignore.case = TRUE) | grepl("eqypt", exp_2016$locationcheck, ignore.case = TRUE), 1, 0)
mean(exp_2016$awarelocation)

log_print("Study 1 Robustness Check, aware of location of videos")
log_print(mean(exp_2016$awarelocation))
###15 respondents have both hk and china in their answers

##subset of respondents unaware of the protest location

unawareloc <-exp_2016[exp_2016$awarelocation ==0,]
nrow(unawareloc)
##robustness check
head(unawareloc)
# RC: bivariate regression
rcH1reg.sanc <- lm(support_sanctions~violence, data=unawareloc)
summary(rcH1reg.sanc)
rcH1reg.don <- lm(money~violence, data=unawareloc)
summary(rcH1reg.don)
rcH1reg.behav <- lm(pol_behavior~violence, data=unawareloc)
summary(rcH1reg.behav)
stargazer(rcH1reg.sanc, rcH1reg.don, rcH1reg.behav, single.row=TRUE, title="Robustness Check for Study 1: Effects of Violent Police Response")

###plot violence effects for all three DVs

RCdiff.in.means_95 <- t.test(unawareloc[unawareloc$violence==1,]$support_sanctions, unawareloc[unawareloc$violence==0,]$support_sanctions,conf.level = .95)
RCdiff.in.means_90 <- t.test(unawareloc[unawareloc$violence==1,]$support_sanctions, unawareloc[unawareloc$violence==0,]$support_sanctions,conf.level = .90)

RCeffect1 <- RCdiff.in.means_95$estimate[1]-RCdiff.in.means_95$estimate[2]

RCdiff.in.means2_95 <- t.test(unawareloc[unawareloc$violence==1,]$money, unawareloc[unawareloc$violence==0,]$money, conf.level = .95)
RCdiff.in.means2_90 <- t.test(unawareloc[unawareloc$violence==1,]$money, unawareloc[unawareloc$violence==0,]$money, conf.level = .90)

RCeffect2 <- RCdiff.in.means2_95$estimate[1]-RCdiff.in.means2_95$estimate[2]

RCdiff.in.means3_95 <- t.test(unawareloc[unawareloc$violence==1,]$pol_behavior, unawareloc[unawareloc$violence==0,]$pol_behavior, conf.level=.95)
RCdiff.in.means3_90 <- t.test(unawareloc[unawareloc$violence==1,]$pol_behavior, unawareloc[unawareloc$violence==0,]$pol_behavior, conf.level=.90)

RCeffect3 <- RCdiff.in.means3_95$estimate[1]-RCdiff.in.means3_95$estimate[2]

##PLOT
dev.off()
par(mfrow=c(3,1))
plot(RCeffect1,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(RCdiff.in.means_95$conf.int[1],3,RCdiff.in.means_95$conf.int[2],3,lwd=2)
segments(RCdiff.in.means_90$conf.int[1],3,RCdiff.in.means_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("RC:Effect of Violence on Support for Sanctions",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(RCeffect2,3,cex=3,pch=20,xlim=c(-15,15),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(RCdiff.in.means2_95$conf.int[1],3,RCdiff.in.means2_95$conf.int[2],3,lwd=2)
segments(RCdiff.in.means2_90$conf.int[1],3,RCdiff.in.means2_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("RC:Effect of Violence on Money Donations",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-15,15,by=5)),labels=c("-15", "-10", "-5", "0", "5", "10", "15"))
mtext("Dollar Amount Change",side=1,cex=1.1,line=2.3)

plot(RCeffect3,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(RCdiff.in.means3_95$conf.int[1],3,RCdiff.in.means3_95$conf.int[2],3,lwd=2)
segments(RCdiff.in.means3_90$conf.int[1],3,RCdiff.in.means3_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("RC:Effect of Violence on Political Behavior",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

## PLOT FOR APPENDIX FIGURE 2
# comparison of main table and robustness check
dev.off()
par(mfrow=c(3,2))
#par(oma=c(0,0,2,2)) # (if following plots don't work, remove/add this line)
plot(effect1,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_95$conf.int[1],3,diff.in.means_95$conf.int[2],3,lwd=2)
segments(diff.in.means_90$conf.int[1],3,diff.in.means_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Support for Sanctions \n all data, N = 1008",side=3,cex=1.1,line=1)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(RCeffect1,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(RCdiff.in.means_95$conf.int[1],3,RCdiff.in.means_95$conf.int[2],3,lwd=2)
segments(RCdiff.in.means_90$conf.int[1],3,RCdiff.in.means_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Support for Sanctions \n no location recognition, N = 736",side=3,cex=1.1,line=1)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect2,3,cex=3,pch=20,xlim=c(-11,11),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(diff.in.means2_95$conf.int[1],3,diff.in.means2_95$conf.int[2],3,lwd=2)
segments(diff.in.means2_90$conf.int[1],3,diff.in.means2_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Money Donations \n all data, N = 1008",side=3,cex=1.1,line=1)
axis(1,at=c(seq(-15,15,by=5)),labels=c("-15", "-10", "-15", "0", "5", "10", "15"))
mtext("Dollar Amount Change",side=1,cex=1.1,line=2.3)

plot(RCeffect2,3,cex=3,pch=20,xlim=c(-15,15),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(RCdiff.in.means2_95$conf.int[1],3,RCdiff.in.means2_95$conf.int[2],3,lwd=2)
segments(RCdiff.in.means2_90$conf.int[1],3,RCdiff.in.means2_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Money Donations \n no location recognition, N = 736",side=3,cex=1.1,line=1)
axis(1,at=c(seq(-15,15,by=5)),labels=c("-15", "-10", "-5", "0", "5", "10", "15"))
mtext("Dollar Amount Change",side=1,cex=1.1,line=2.3)

plot(effect3,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(diff.in.means3_95$conf.int[1],3,diff.in.means3_95$conf.int[2],3,lwd=2)
segments(diff.in.means3_90$conf.int[1],3,diff.in.means3_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Political Behavior  \n all data, N = 1008",side=3,cex=1.1,line=1)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(RCeffect3,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(RCdiff.in.means3_95$conf.int[1],3,RCdiff.in.means3_95$conf.int[2],3,lwd=2)
segments(RCdiff.in.means3_90$conf.int[1],3,RCdiff.in.means3_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Political Behavior \n no location recognition, N = 736",side=3,cex=1.1,line=1)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

##################################################################
##################################################################
###################### STUDY 2 - 2018 ############################
##################################################################
##################################################################

## EXPLANATION: LINES 899 TO 1014 LOAD THE DATA FOR STUDY 2 IN THE PAPER. ##


# setwd() -- SET TO APPROPRIATE FOLDER

exp <- read.csv("FGK_study2.csv")

##remove the first row which lists question name
exp <- exp[-1,]

##drop respondents who did not consent
exp <- exp[exp$presurvey==1,]

exp <- exp[exp$Q165 == 1,] #subset to those who completed
nrow(exp)
## DEMOGRAPHICS

## AGE
exp$age <- 2018 - as.numeric(as.character(exp$birthY))
## GENDER
exp$male <- ifelse(exp$gender==1, 1, 0)
## ETHNICITY
exp$ethnicity_1 <- ifelse(exp$ethnicity_1 == 1, 1, 0)
exp$ethnicity_2 <- ifelse(exp$ethnicity_2 == 1, 1, 0)
exp$ethnicity_3 <- ifelse(exp$ethnicity_3 == 1, 1, 0)
exp$ethnicity_4 <- ifelse(exp$ethnicity_4 == 1, 1, 0)
exp$ethnicity_5 <- ifelse(exp$ethnicity_5 == 1, 1, 0)
exp$ethnicity_6 <- ifelse(exp$ethnicity_6 == 1, 1, 0)
exp$ethnicity_7 <- ifelse(exp$ethnicity_7 == 1, 1, 0)
exp$ethnicity_sum <- as.numeric(as.character(exp$ethnicity_1)) + as.numeric(as.character(exp$ethnicity_2)) + as.numeric(as.character(exp$ethnicity_3)) + as.numeric(as.character(exp$ethnicity_4)) + as.numeric(as.character(exp$ethnicity_5)) + as.numeric(as.character(exp$ethnicity_6)) + as.numeric(as.character(exp$ethnicity_7))
exp$ethnicity <- ifelse(exp$ethnicity_sum > 1, "multiracial", NA)
exp$ethnicity <- ifelse(exp$ethnicity_sum == 1 & exp$ethnicity_1 == 1, "black", exp$ethnicity)
exp$ethnicity <- ifelse(exp$ethnicity_sum == 1 & exp$ethnicity_2 == 1, "white", exp$ethnicity)
exp$ethnicity <- ifelse(exp$ethnicity_sum == 1 & exp$ethnicity_3 == 1, "hispanic", exp$ethnicity)
exp$ethnicity <- ifelse(exp$ethnicity_sum == 1 & exp$ethnicity_4 == 1, "asian", exp$ethnicity)
exp$ethnicity <- ifelse(exp$ethnicity_sum == 1 & exp$ethnicity_5 == 1, "arab", exp$ethnicity)
exp$ethnicity <- ifelse(exp$ethnicity_sum == 1 & exp$ethnicity_6 == 1, "native_am", exp$ethnicity)
exp$ethnicity <- ifelse(exp$ethnicity_sum == 1 & exp$ethnicity_7 == 1, "other", exp$ethnicity)
exp$white <- ifelse(exp$ethnicity == "white", 1, 0)
## PROTESTER
exp$protester <- ifelse(exp$protest_1 == 1 | exp$protest_2 == 1 | exp$protest_3 == 1 | exp$protest_4 == 1 | exp$protest_5 == 1 | exp$protest_6 == 1 | exp$protest_7 == 1 | exp$protest_8 == 1, 1, 0)
exp$trump_era_protester <- ifelse(exp$protest_1 == 1, 1, 0)
## INCOME
exp$mid_income <- ifelse(as.numeric(as.character(exp$income)) >= 5, 1, 0)
exp$mid_income
## PARTY
exp$republican <- ifelse(exp$Q204==2 | exp$partyID==2 | exp$leaning==1 | exp$Q205==1, 1, 0)
exp$democrat <- ifelse(exp$Q204==1 | exp$partyID==1 | exp$leaning==2 | exp$Q205==2, 1, 0)

## VOTE AND VOTE CHOICE
exp$trump <- ifelse(exp$votechoice==2 | exp$Q207==2 , 1, 0)
exp$vote_presidential <- ifelse(exp$Q206==3 | exp$vote==3, 1, 0)
## EDUCATION
exp$college_degree <- ifelse(as.numeric(as.character(exp$edu))>=5, 1, 0)

# other controls
## FEELING THERMOMETER
exp$Latin_thermometer <- as.numeric(as.character(exp$thermometer_5))
## CONFIDENCE IN INSTITUTIONS
exp$USmedia_confidence <- as.numeric(as.character(exp$confidence_1))
exp$USgovt_confidence <- as.numeric(as.character(exp$confidence_3))
exp$USpolice_confidence <- as.numeric(as.character(exp$confidence_4))

##attention to int'l news
exp$attention_news <- car::recode(exp$news, "'5'='0'; '4'='1'; '3'='2'; '2'='3'; '1'='4'; 'NA'='NA'")
##news from social media
exp$social_media_news <- car::recode(exp$news_source_1, "'2'='1'; '3'='1'; '1'=0; 'NA'='NA'" )
exp$social_media_news <- as.numeric(as.character(exp$social_media_news))
##news from tv
exp$tv_news <- car::recode(exp$news_source_6, "'2'='1'; '3'='1'; '1'=0; 'NA'='NA'")
exp$tv_news <- as.numeric(as.character(exp$tv_news))

###TREATMENTS
exp$phonenonviolent <- ifelse(exp$DO.BL.video=="timer|phonenonviolent", 1, 0)
sum(exp$phonenonviolent) #460
exp$phoneviolent <- ifelse(exp$DO.BL.video=="timer|phoneviolent", 1, 0)
sum(exp$phoneviolent) #425
exp$newsviolent <- ifelse(exp$DO.BL.video=="timer|newsviolent", 1, 0)
sum(exp$newsviolent) # 455
exp$newsnonviolent <- ifelse(exp$DO.BL.video=="timer|newsnonviolent", 1, 0)
sum(exp$newsnonviolent) #408
exp$mobile <- ifelse(exp$phonenonviolent==1 | exp$phoneviolent==1, 1, 0)
exp$violence <- ifelse(exp$newsviolent==1 | exp$phoneviolent==1, 1, 0)

###OUTCOMES
##sanctions
exp$sanctions.ordinal <- ifelse(exp$sanctions.preference=="1",5, ifelse(exp$sanctions.preference=="4",4, ifelse(exp$sanctions.preference=="5",3,ifelse(exp$sanctions.preference=="6",2,1))))
exp$sanctions.binary <- ifelse(exp$sanctions.preference=="1"|exp$sanctions.preference=="4", 1, 0)
exp$sanctions.unsure <- ifelse(exp$sanctions.preference == "5", 1, 0)
##learnmore
exp$learnmore <- ifelse(exp$learnmore=="1", 1, 0)

#how many respondents want to learn more? -- footnote 10 in main paper
mean(exp[exp$violence ==1 & exp$mobile ==1,]$learnmore)
mean(exp[exp$violence ==0 & exp$mobile ==1,]$learnmore)
mean(exp[exp$violence ==0 & exp$mobile ==0,]$learnmore)
mean(exp[exp$violence ==1 & exp$mobile ==0,]$learnmore)

log_print("Study 2, how many respondents want to learn more in each treatment condition?")

log_print(mean(exp[exp$violence ==1 & exp$mobile ==1,]$learnmore))
log_print(mean(exp[exp$violence ==0 & exp$mobile ==1,]$learnmore))
log_print(mean(exp[exp$violence ==0 & exp$mobile ==0,]$learnmore))
log_print(mean(exp[exp$violence ==1 & exp$mobile ==0,]$learnmore))

##blame of protesters
exp$blame_2 <- as.numeric(as.character(exp$blame_2))
##blame of police
exp$blame_1 <- as.numeric(as.character(exp$blame_1))
##blame of government
exp$blame_3 <- as.numeric(as.character(exp$blame_3))

exp$blameprotestor <- exp$blame_2
exp$blamepolice <- exp$blame_1
exp$blamegov<-exp$blame_3

##trustworthy
exp$trustworthy.ordinal <- ifelse(exp$trustvideo=="1",5, ifelse(exp$trustvideo=="2",4, ifelse(exp$trustvideo=="4",3,ifelse(exp$trustvideo=="6",2, ifelse(exp$trustvideo=="7",1,0)))))
exp$trustworthy.binary <- ifelse(exp$trustvideo=="1"|exp$trustvideo=="2"|exp$trustvideo=="4", 1, 0)

#########################Analysis##################################

## EXPLANATION: LINES 1038 TO 2465 GIVE RESULTS FOR STUDY 2 IN THE PAPER ##
## RUN STUDY 2 CODE IN FULL (LINES 899 TO 2465) TO MAKE SURE EVERYTHING WORKS AS CODE OFTEN RELIES ON CODE BEFORE IT ##

## EXPLANATION: THE MEANS IN THE PAPER (WRITTEN IN TEXT ON PAGE 12) CAN BE FOUND IN LINES 1048 to 1054 ##

## EXPLANATION: LINES 1092 TO 1112 GIVES PAPER FIGURE 4 PLOT ### 

## EXPLANATION: LINES 1284 TO 1312 GIVES PAPER FIGURE 5 PLOT ### 

## EXPLANATION: LINES 1355 TO 1382 GIVES PAPER FIGURE 6 PLOT ### 

## EXPLANATION: THE ~30% ESTIMATE IN THE PAPER (ON PAGE 13, WRITTEN IN TEXT OF FOOTNOTE 10) CAN BE FOUND IN LINES 988 to 999 ##

## EXPLANATION: THE TRUSTWORTHINESS ESTIMATES IN THE PAPER (ON PAGE 14) CAN BE FOUND IN LINES 1384 to 1392 ##

## EXPLANATION: LINES 2918 TO 2923 GIVE APPENDIX TABLE 8 (STUDY 2 REGRESSION TABLE) ##

## EXPLANATION: LINES 2925 TO 2930 GIVE APPENDIX TABLE 9 (STUDY 2 REGRESSION TABLE) ##


## T-TESTS AND REGRESSIONS

## Regressions WITHOUT Interaction term, violence

## Support for Sanctions
##check difference in means for effect of violence 

diff.in.means_sanc_violence_95 <- t.test(exp[exp$violence==1,]$sanctions.binary, exp[exp$violence==0,]$sanctions.binary)
diff.in.means_sanc_violence_90 <- t.test(exp[exp$violence==1,]$sanctions.binary, exp[exp$violence==0,]$sanctions.binary, conf.level = 0.9)

# mean for nonviolent videos

diff.in.means_sanc_violence_95$estimate[2]

# mean for violent videos

diff.in.means_sanc_violence_95$estimate[1] 

effect_sanc_violence <- diff.in.means_sanc_violence_95$estimate[1]-diff.in.means_sanc_violence_95$estimate[2]

log_print("Study 2, support for sanctions, difference in means")
log_print(diff.in.means_sanc_violence_95)

# bivariate regression with violence as IV - violence is significant, p = 0.0036

reg.binary_sanc.violence.bivar <- lm(sanctions.binary ~ violence, data=exp)
summary(reg.binary_sanc.violence.bivar)

##regression with all covariates - violence is significant, p = 0.002297

reg.binary_sanc.violence.covar <- lm(sanctions.binary ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.binary_sanc.violence.covar)


## Learning More
##check difference in means for effect of violence 

diff.in.means_learn_violence_95 <- t.test(exp[exp$violence==1,]$learnmore, exp[exp$violence==0,]$learnmore)
diff.in.means_learn_violence_90 <- t.test(exp[exp$violence==1,]$learnmore, exp[exp$violence==0,]$learnmore, conf.level = 0.9)

effect_learn_violence <- diff.in.means_learn_violence_95$estimate[1]-diff.in.means_learn_violence_95$estimate[2]

# bivariate regression with violence as IV - violence is not significant

reg.learn.violence.bivar <- lm(learnmore ~ violence, data=exp)
summary(reg.learn.violence.bivar)

##regression with all covariates - violence is not significant

reg.learn.violence.covar <- lm(learnmore ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.learn.violence.covar)

###plot violence effects on sanctions support and learning more in regression without interaction

## PLOT FOR FIGURE 4
## Difference in means plots

#pdf("study2_violence.pdf")
par(mfrow=c(2,1))
plot(effect_sanc_violence,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_sanc_violence_95$conf.int[1],3,diff.in.means_sanc_violence_95$conf.int[2],3,lwd=2)
segments(diff.in.means_sanc_violence_90$conf.int[1],3,diff.in.means_sanc_violence_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Support for Sanctions",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect_learn_violence,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_learn_violence_95$conf.int[1],3,diff.in.means_learn_violence_95$conf.int[2],3,lwd=2)
segments(diff.in.means_learn_violence_90$conf.int[1],3,diff.in.means_learn_violence_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Interest in Learning More",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
#dev.off()


## Regressions WITH Interaction term

## Support for Sanctions
##check difference in means for effect of violence 

# trivariate regression - violence is significant, p = 0.0139, mobile is not significant, neither is interaction

reg.binary_sanc.violencemobile.inter <- lm(sanctions.binary ~  mobile + violence + mobile*violence, data=exp)
summary(reg.binary_sanc.violencemobile.inter)

# reversing coding of video source (so that mobile = 0 and tv = 1) shows that violence would also be significant in that equation
exp$TV <- ifelse(exp$mobile == 1, 0, 1)

# reversed coding of source, trivariate regression - violence is still significant, p = 0.0871, mobile is not significant, neither is interaction
reg.binary_sanc.violencemobile.inter_rev <- lm(sanctions.binary ~  TV + violence + TV*violence, data=exp)
summary(reg.binary_sanc.violencemobile.inter_rev)

#calculating marginal effect of mobile when violence = 1 for sanctions support (binary coding of sanctions support)

vcov(reg.binary_sanc.violencemobile.inter)[4,2]
# var for mobile is 0.0011511991, vcov(reg.binary_sanc.violencemobile.inter)[2,2]
# var for interaction term is 0.0022839411, vcov(reg.binary_sanc.violencemobile.inter)[4,4]
# cov for mobile and interaction is -0.0011511991, vcov(reg.binary_sanc.violencemobile.inter)[4,2]
margeffect.bin.sanc <- summary(reg.binary_sanc.violencemobile.inter)$coef[2] + summary(reg.binary_sanc.violencemobile.inter)$coef[4]
var_margeffect.bin.sanc <- vcov(reg.binary_sanc.violencemobile.inter)[2,2] + vcov(reg.binary_sanc.violencemobile.inter)[4,4] + (2 * vcov(reg.binary_sanc.violencemobile.inter)[4,2])
se_margeffect.bin.sanc <- sqrt(var_margeffect.bin.sanc)

ci_upper.bin.sanc <- margeffect.bin.sanc + (1.96*se_margeffect.bin.sanc)
ci_lower.bin.sanc <- margeffect.bin.sanc - (1.96*se_margeffect.bin.sanc)

ci_upper.bin.sanc90 <- margeffect.bin.sanc + (1.645*se_margeffect.bin.sanc)
ci_lower.bin.sanc90 <- margeffect.bin.sanc - (1.645*se_margeffect.bin.sanc)


## marginal effect is not significant, as confidence interval crosses zero


#calculating marginal effect of violence when mobile = 1 for sanctions support (binary coding of sanctions support)

vcov(reg.binary_sanc.violencemobile.inter)[4,3]
# var for violence is 0.001157145, vcov(reg.binary_sanc.violencemobile.inter)[3,3]
# var for interaction term is 0.0022839411, vcov(reg.binary_sanc.violencemobile.inter)[4,4]
# cov for mobile and interaction is -0.001157145, vcov(reg.binary_sanc.violencemobile.inter)[4,3]

margeffect.bin.sanc.vio <- summary(reg.binary_sanc.violencemobile.inter)$coef[3] + summary(reg.binary_sanc.violencemobile.inter)$coef[4]
var_margeffect.bin.sanc.vio <- vcov(reg.binary_sanc.violencemobile.inter)[3,3] + vcov(reg.binary_sanc.violencemobile.inter)[4,4] + (2 * vcov(reg.binary_sanc.violencemobile.inter)[4,3])
se_margeffect.bin.sanc.vio <- sqrt(var_margeffect.bin.sanc.vio)

ci_upper.bin.sanc.vio <- margeffect.bin.sanc.vio + (1.96*se_margeffect.bin.sanc.vio)
ci_lower.bin.sanc.vio <- margeffect.bin.sanc.vio - (1.96*se_margeffect.bin.sanc.vio)

ci_upper.bin.sanc90.vio <- margeffect.bin.sanc.vio + (1.645*se_margeffect.bin.sanc.vio)
ci_lower.bin.sanc90.vio <- margeffect.bin.sanc.vio - (1.645*se_margeffect.bin.sanc.vio)


## marginal effect is significant for 90% CI, not for 95% CI


# # 
# ##regression with interaction and all covariates - violence is significant, p = 0.0009502

reg.binary_sanc.violencemobile.covar.inter <- lm(sanctions.binary ~ mobile + violence + mobile*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.binary_sanc.violencemobile.covar.inter)

# # reversing again, just to check, in the regression for all covariates, violence is still significant, p = 0.084985
# reg.binary_sanc.violencemobile.covar.inter_rev <- lm(sanctions.binary ~ TV + violence + TV*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
# summary(reg.binary_sanc.violencemobile.covar.inter_rev)


# ## Learning More

# trivariate regression - violence is not significant, mobile is not significant, neither is interaction

reg.learn.violencemobile.inter <- lm(learnmore ~ mobile + violence +  mobile*violence, data=exp)
summary(reg.learn.violencemobile.inter)

#calculating marginal effect of mobile when violence = 1 for learning more

vcov(reg.learn.violencemobile.inter)
# var for mobile is 0.0009690888, vcov(reg.learn.violencemobile.inter)[2,2]
# var for interaction term is 0.0019226403, vcov(reg.learn.violencemobile.inter)[4,4]
# cov for mobile and interaction is -0.0009690888, vcov(reg.learn.violencemobile.inter)[4,2]


margeffect.bin.learnmore <- summary(reg.learn.violencemobile.inter)$coef[2] + summary(reg.learn.violencemobile.inter)$coef[4]
var_margeffect.bin.learnmore <- vcov(reg.learn.violencemobile.inter)[2,2] + vcov(reg.learn.violencemobile.inter)[4,4] + (2 * vcov(reg.learn.violencemobile.inter)[4,2])
se_margeffect.bin.learnmore <- sqrt(var_margeffect.bin.learnmore)

ci_upper.bin.learnmore <- margeffect.bin.learnmore + (1.96*se_margeffect.bin.learnmore)
ci_lower.bin.learnmore <- margeffect.bin.learnmore - (1.96*se_margeffect.bin.learnmore)

ci_upper.bin.learnmore90 <- margeffect.bin.learnmore + (1.645*se_margeffect.bin.learnmore)
ci_lower.bin.learnmore90 <- margeffect.bin.learnmore - (1.645*se_margeffect.bin.learnmore)


## marginal effect is not significant, as confidence interval crosses zero for both CIs



#calculating marginal effect of violence when mobile = 1 for learning more

vcov(reg.learn.violencemobile.inter)
# var for violence is 0.0009740944, vcov(reg.learn.violencemobile.inter)[3,3]
# var for interaction term is 0.0019226403, vcov(reg.learn.violencemobile.inter)[4,4]
# cov for mobile and interaction is -0.0009740944, vcov(reg.learn.violencemobile.inter)[4,3]


margeffect.bin.learnmore.vio <- summary(reg.learn.violencemobile.inter)$coef[3] + summary(reg.learn.violencemobile.inter)$coef[4]
var_margeffect.bin.learnmore.vio <- vcov(reg.learn.violencemobile.inter)[3,3] + vcov(reg.learn.violencemobile.inter)[4,4] + (2 * vcov(reg.learn.violencemobile.inter)[4,3])
se_margeffect.bin.learnmore.vio <- sqrt(var_margeffect.bin.learnmore.vio)

ci_upper.bin.learnmore.vio <- margeffect.bin.learnmore.vio + (1.96*se_margeffect.bin.learnmore.vio)
ci_lower.bin.learnmore.vio <- margeffect.bin.learnmore.vio - (1.96*se_margeffect.bin.learnmore.vio)

ci_upper.bin.learnmore90.vio <- margeffect.bin.learnmore.vio + (1.645*se_margeffect.bin.learnmore.vio)
ci_lower.bin.learnmore90.vio <- margeffect.bin.learnmore.vio - (1.645*se_margeffect.bin.learnmore.vio)


## marginal effect is not significant, as confidence interval crosses zero for both CIs

##regression with interaction and all covariates - violence is not significant, mobile is not significant, neither is interaction

reg.learn.violencemobile.covar.inter <- lm(learnmore ~ mobile + violence + mobile*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.learn.violencemobile.covar.inter)


###plot interaction effects on support for sanctions

## Difference in means plots

# sanctions support

#binary_sanc.inter.mobile <- summary(reg.binary_sanc.violencemobile.inter)$coefficients[2]
binary_sanc.inter.violence <- summary(reg.binary_sanc.violencemobile.inter)$coefficients[3]
binary_sanc.inter.margeffect.vio <- margeffect.bin.sanc.vio
binary_sanc.inter.margeffect <- margeffect.bin.sanc

#binary_sanc.inter.mobile.err <- summary(reg.binary_sanc.violencemobile.inter)$coefficients[,2][2]
#li_binary_sanc.inter.mobile <- binary_sanc.inter.mobile-1.96 *binary_sanc.inter.mobile.err
#ui_binary_sanc.inter.mobile <- binary_sanc.inter.mobile+1.96 *binary_sanc.inter.mobile.err

binary_sanc.inter.violence.err <- summary(reg.binary_sanc.violencemobile.inter)$coefficients[,2][3]
li_binary_sanc.inter.violence <- binary_sanc.inter.violence-1.96 *binary_sanc.inter.violence.err
ui_binary_sanc.inter.violence <- binary_sanc.inter.violence+1.96 *binary_sanc.inter.violence.err
li_binary_sanc.inter.violence90 <- binary_sanc.inter.violence-1.645 *binary_sanc.inter.violence.err
ui_binary_sanc.inter.violence90 <- binary_sanc.inter.violence+1.645 *binary_sanc.inter.violence.err


li_binary_sanc.inter.margeffect.vio <- ci_lower.bin.sanc.vio
ui_binary_sanc.inter.margeffect.vio <- ci_upper.bin.sanc.vio
li_binary_sanc.inter.margeffect90.vio <- ci_lower.bin.sanc90.vio
ui_binary_sanc.inter.margeffect90.vio <- ci_upper.bin.sanc90.vio


li_binary_sanc.inter.margeffect <- ci_lower.bin.sanc
ui_binary_sanc.inter.margeffect <- ci_upper.bin.sanc
li_binary_sanc.inter.margeffect90 <- ci_lower.bin.sanc90
ui_binary_sanc.inter.margeffect90 <- ci_upper.bin.sanc90


binary_sanc.inter_coef <- c(binary_sanc.inter.violence, binary_sanc.inter.margeffect.vio, binary_sanc.inter.margeffect)
binary_sanc.inter_li <- c(li_binary_sanc.inter.violence, li_binary_sanc.inter.margeffect.vio, li_binary_sanc.inter.margeffect)
binary_sanc.inter_ui <- c(ui_binary_sanc.inter.violence, ui_binary_sanc.inter.margeffect.vio, ui_binary_sanc.inter.margeffect)
binary_sanc.inter_li90 <- c(li_binary_sanc.inter.violence90, li_binary_sanc.inter.margeffect90.vio, li_binary_sanc.inter.margeffect90)
binary_sanc.inter_ui90 <- c(ui_binary_sanc.inter.violence90, ui_binary_sanc.inter.margeffect90.vio, ui_binary_sanc.inter.margeffect90)

treatment <- c("violence, m=0", "violence, m=1", "mobile, v=1")
df_binary_sanc.inter <- as.data.frame(cbind(treatment, binary_sanc.inter_coef, binary_sanc.inter_li90, binary_sanc.inter_ui90, binary_sanc.inter_li, binary_sanc.inter_ui))

# PLOT FOR FIGURE 5
#pdf("study2_sanctions_new_w90.pdf")
par(mfrow=c(3,1))

plot(binary_sanc.inter.violence,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_sanc.inter.violence,3,ui_binary_sanc.inter.violence,3,lwd=2)
segments(li_binary_sanc.inter.violence90,3,ui_binary_sanc.inter.violence90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.9,line=2.3)

plot(binary_sanc.inter.margeffect.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(li_binary_sanc.inter.margeffect.vio,3,ui_binary_sanc.inter.margeffect.vio,3,lwd=2)
segments(li_binary_sanc.inter.margeffect90.vio,3,ui_binary_sanc.inter.margeffect90.vio,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.9,line=2.3)


plot(binary_sanc.inter.margeffect,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_sanc.inter.margeffect ,3,ui_binary_sanc.inter.margeffect ,3,lwd=2)
segments(li_binary_sanc.inter.margeffect90 ,3,ui_binary_sanc.inter.margeffect90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.9,line=2.3)
#dev.off()


# learning more

#learn.inter.mobile <- summary(reg.learn.violencemobile.inter)$coefficients[2]
learn.inter.violence <- summary(reg.learn.violencemobile.inter)$coefficients[3]
learn.inter.margeffect.vio <- margeffect.bin.learnmore.vio
learn.inter.margeffect <- margeffect.bin.learnmore

#learn.inter.mobile.err <- summary(reg.learn.violencemobile.inter)$coefficients[,2][2]
#li_learn.inter.mobile <- learn.inter.mobile-1.96 *learn.inter.mobile.err
#ui_learn.inter.mobile <- learn.inter.mobile+1.96 *learn.inter.mobile.err

learn.inter.violence.err <- summary(reg.learn.violencemobile.inter)$coefficients[,2][3]
li_learn.inter.violence <- learn.inter.violence-1.96 *learn.inter.violence.err
ui_learn.inter.violence<- learn.inter.violence+1.96 *learn.inter.violence.err
li_learn.inter.violence90 <- learn.inter.violence-1.645 *learn.inter.violence.err
ui_learn.inter.violence90 <- learn.inter.violence+1.645 *learn.inter.violence.err


li_learn.inter.margeffect.vio <- ci_lower.bin.learnmore.vio
ui_learn.inter.margeffect.vio <- ci_upper.bin.learnmore.vio
li_learn.inter.margeffect90.vio <- ci_lower.bin.learnmore90.vio
ui_learn.inter.margeffect90.vio <- ci_upper.bin.learnmore90.vio



li_learn.inter.margeffect <- ci_lower.bin.learnmore
ui_learn.inter.margeffect <- ci_upper.bin.learnmore
li_learn.inter.margeffect90 <- ci_lower.bin.learnmore90
ui_learn.inter.margeffect90 <- ci_upper.bin.learnmore90


learn.inter_coef <- c(learn.inter.violence, learn.inter.margeffect.vio, learn.inter.margeffect)
learn.inter_li <- c(li_learn.inter.violence, li_learn.inter.margeffect.vio, li_learn.inter.margeffect)
learn.inter_ui <- c(ui_learn.inter.violence, ui_learn.inter.margeffect.vio, ui_learn.inter.margeffect)
learn.inter_li90 <- c(li_learn.inter.violence90, li_learn.inter.margeffect90.vio, li_learn.inter.margeffect90)
learn.inter_ui90 <- c(ui_learn.inter.violence90, ui_learn.inter.margeffect90.vio, ui_learn.inter.margeffect90)

treatment <- c("violence, m=0", "violence, m=1", "mobile, v=1")
df_learn.inter <- as.data.frame(cbind(treatment, learn.inter_coef, learn.inter_li90, learn.inter_ui90, learn.inter_li, learn.inter_ui))

# PLOT FOR FIGURE 6
#pdf("study2_learn_more_new_w90.pdf")
par(mfrow=c(3,1))

plot(learn.inter.violence,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_learn.inter.violence,3,ui_learn.inter.violence,3,lwd=2)
segments(li_learn.inter.violence90,3,ui_learn.inter.violence90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.9,line=2.3)

plot(learn.inter.margeffect.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(li_learn.inter.margeffect.vio,3,ui_learn.inter.margeffect.vio,3,lwd=2)
segments(li_learn.inter.margeffect90.vio,3,ui_learn.inter.margeffect90.vio,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.9,line=2.3)

plot(learn.inter.margeffect,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_learn.inter.margeffect ,3,ui_learn.inter.margeffect ,3,lwd=2)
segments(li_learn.inter.margeffect90 ,3,ui_learn.inter.margeffect90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.9,line=2.3)
#dev.off()

# trustworthiness of mobile/tv news clips - mentioned on page 14

mean(exp[exp$mobile ==1,]$trustworthy.ordinal)
mean(exp[exp$mobile ==0,]$trustworthy.ordinal)

log_print("Study 2, trustworthiness of mobile clips")
log_print(mean(exp[exp$mobile ==1,]$trustworthy.ordinal))
log_print("Study 2, trustworthiness of TV clips")
log_print(mean(exp[exp$mobile ==0,]$trustworthy.ordinal))
#######################SUPPLEMENTAL#####
##verify that results using ordinal outcome variables are consistent with findings - statement p. 6 of paper (results not illustrated in paper or appendix)

## Regressions WITH Interaction term, ordinal sanctions support

## Support for Sanctions
##check difference in means for effect of violence 

# trivariate regression - violence is still significant, p = 0.0114, mobile is still not significant, neither is interaction (still)

reg.ord_sanc.violencemobile.inter <- lm(sanctions.ordinal ~  mobile + violence + mobile*violence, data=exp)
summary(reg.ord_sanc.violencemobile.inter)

# reversing coding of video source (so that mobile = 0 and tv = 1) shows that violence would also be significant in that equation

# reversed coding of source, trivariate regression - violence is still significant, p = 0.0363, mobile is not significant, neither is interaction
reg.ord_sanc.violencemobile.inter_rev <- lm(sanctions.ordinal ~  TV + violence + TV*violence, data=exp)
summary(reg.ord_sanc.violencemobile.inter_rev)

#calculating marginal effect of mobile when violence = 1 for sanctions support (ordinal coding of sanctions support)

margeffect.ord.sanc <- summary(reg.ord_sanc.violencemobile.inter)$coef[2] + summary(reg.ord_sanc.violencemobile.inter)$coef[4]
var_margeffect.ord.sanc <- vcov(reg.ord_sanc.violencemobile.inter)[2,2] + vcov(reg.ord_sanc.violencemobile.inter)[4,4] + (2 * vcov(reg.ord_sanc.violencemobile.inter)[4,2])
se_margeffect.ord.sanc <- sqrt(var_margeffect.ord.sanc)

ci_upper.ord.sanc <- margeffect.ord.sanc + (1.96*se_margeffect.ord.sanc)
ci_lower.ord.sanc <- margeffect.ord.sanc - (1.96*se_margeffect.ord.sanc)

ci_upper.ord.sanc90 <- margeffect.ord.sanc + (1.645*se_margeffect.ord.sanc)
ci_lower.ord.sanc90 <- margeffect.ord.sanc - (1.645*se_margeffect.ord.sanc)

#calculating marginal effect of violence when mobile = 1 for sanctions support (ordinal coding of sanctions support)

margeffect.ord.sanc.vio <- summary(reg.ord_sanc.violencemobile.inter)$coef[3] + summary(reg.ord_sanc.violencemobile.inter)$coef[4]
var_margeffect.ord.sanc.vio <- vcov(reg.ord_sanc.violencemobile.inter)[3,3] + vcov(reg.ord_sanc.violencemobile.inter)[4,4] + (2 * vcov(reg.ord_sanc.violencemobile.inter)[4,3])
se_margeffect.ord.sanc.vio <- sqrt(var_margeffect.ord.sanc.vio)

ci_upper.ord.sanc.vio <- margeffect.ord.sanc.vio + (1.96*se_margeffect.ord.sanc.vio)
ci_lower.ord.sanc.vio <- margeffect.ord.sanc.vio - (1.96*se_margeffect.ord.sanc.vio)

ci_upper.ord.sanc90.vio <- margeffect.ord.sanc.vio + (1.645*se_margeffect.ord.sanc.vio)
ci_lower.ord.sanc90.vio <- margeffect.ord.sanc.vio - (1.645*se_margeffect.ord.sanc.vio)


## marginal effect is still significant, now at both 90% and 95% CIs

##regression with all covariates, ordinal sanctions support coding - violence is still significant, p = 0.00612

reg.ord_sanc.violencemobile.covar.inter <- lm(sanctions.ordinal ~ mobile + violence + mobile*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.ord_sanc.violencemobile.covar.inter)

# reversing again, just to check, in the regression for all covariates, violence is still significant, p = 0.03444
reg.ord_sanc.violencemobile.covar.inter_rev <- lm(sanctions.ordinal ~ TV + violence + TV*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.ord_sanc.violencemobile.covar.inter_rev)

###plot interaction effects on support for sanctions

## Difference in means plots

# sanctions support

#ord_sanc.inter.mobile <- summary(reg.ord_sanc.violencemobile.inter)$coefficients[2]
ord_sanc.inter.violence <- summary(reg.ord_sanc.violencemobile.inter)$coefficients[3]
ord_sanc.inter.margeffect.vio <- margeffect.ord.sanc.vio
ord_sanc.inter.margeffect <- margeffect.ord.sanc

#ord_sanc.inter.mobile.err <- summary(reg.ord_sanc.violencemobile.inter)$coefficients[,2][2]
#li_ord_sanc.inter.mobile <- ord_sanc.inter.mobile-1.96 *ord_sanc.inter.mobile.err
#ui_ord_sanc.inter.mobile <- ord_sanc.inter.mobile+1.96 *ord_sanc.inter.mobile.err

ord_sanc.inter.violence.err <- summary(reg.ord_sanc.violencemobile.inter)$coefficients[,2][3]
li_ord_sanc.inter.violence <- ord_sanc.inter.violence-1.96 *ord_sanc.inter.violence.err
ui_ord_sanc.inter.violence <- ord_sanc.inter.violence+1.96 *ord_sanc.inter.violence.err
li_ord_sanc.inter.violence90 <- ord_sanc.inter.violence-1.645 *ord_sanc.inter.violence.err
ui_ord_sanc.inter.violence90 <- ord_sanc.inter.violence+1.645 *ord_sanc.inter.violence.err

li_ord_sanc.inter.margeffect.vio <- ci_lower.ord.sanc.vio
ui_ord_sanc.inter.margeffect.vio <- ci_upper.ord.sanc.vio
li_ord_sanc.inter.margeffect90.vio <- ci_lower.ord.sanc90.vio
ui_ord_sanc.inter.margeffect90.vio <- ci_upper.ord.sanc90.vio


li_ord_sanc.inter.margeffect <- ci_lower.ord.sanc
ui_ord_sanc.inter.margeffect <- ci_upper.ord.sanc
li_ord_sanc.inter.margeffect90 <- ci_lower.ord.sanc90
ui_ord_sanc.inter.margeffect90 <- ci_upper.ord.sanc90

ord_sanc.inter_coef <- c(ord_sanc.inter.violence, ord_sanc.inter.margeffect.vio, ord_sanc.inter.margeffect)
ord_sanc.inter_li <- c(li_ord_sanc.inter.violence, li_ord_sanc.inter.margeffect.vio, li_ord_sanc.inter.margeffect)
ord_sanc.inter_ui <- c(ui_ord_sanc.inter.violence, ui_ord_sanc.inter.margeffect.vio, ui_ord_sanc.inter.margeffect)
ord_sanc.inter_li90 <- c(li_ord_sanc.inter.violence90, li_ord_sanc.inter.margeffect90.vio, li_ord_sanc.inter.margeffect90)
ord_sanc.inter_ui90 <- c(ui_ord_sanc.inter.violence90, ui_ord_sanc.inter.margeffect90.vio, ui_ord_sanc.inter.margeffect90)

treatment <- c("violence, m=0", "violence, m=1", "mobile, v=1")
df_ord_sanc.inter <- as.data.frame(cbind(treatment, ord_sanc.inter_coef, ord_sanc.inter_li90, ord_sanc.inter_ui90, ord_sanc.inter_li, ord_sanc.inter_ui))


par(mfrow=c(3,1))

plot(ord_sanc.inter.violence,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_ord_sanc.inter.violence,3,ui_ord_sanc.inter.violence,3,lwd=2)
segments(li_ord_sanc.inter.violence90,3,ui_ord_sanc.inter.violence90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(ord_sanc.inter.margeffect.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(li_ord_sanc.inter.margeffect.vio,3,ui_ord_sanc.inter.margeffect.vio,3,lwd=2)
segments(li_ord_sanc.inter.margeffect90.vio,3,ui_ord_sanc.inter.margeffect90.vio,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)


plot(ord_sanc.inter.margeffect,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_ord_sanc.inter.margeffect ,3,ui_ord_sanc.inter.margeffect ,3,lwd=2)
segments(li_ord_sanc.inter.margeffect90 ,3,ui_ord_sanc.inter.margeffect90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

# same results, except now effect of violence conditional on M=1 is siginficant at both CIs



###check effects when DVs are ordinal (only applicable for sanctions DV)

## Regressions WITHOUT Interaction term, violence

## Support for Sanctions
##check difference in means for effect of violence (still significant)

diff.in.means_sanc_violence.ord_95 <- t.test(exp[exp$violence==1,]$sanctions.ordinal, exp[exp$violence==0,]$sanctions.ordinal)
diff.in.means_sanc_violence.ord_90 <- t.test(exp[exp$violence==1,]$sanctions.ordinal, exp[exp$violence==0,]$sanctions.ordinal, conf.level = 0.9)

effect_sanc_violence.ord <- diff.in.means_sanc_violence.ord_95$estimate[1]-diff.in.means_sanc_violence.ord_95$estimate[2]

# bivariate regression with violence as IV - violence is still significant, p = 0.00128

reg.ord_sanc.violence.bivar <- lm(sanctions.ordinal ~ violence, data=exp)
summary(reg.ord_sanc.violence.bivar)

##regression with all covariates - violence is still significant, p = 0.000601

reg.ord_sanc.violence.covar <- lm(sanctions.ordinal ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.ord_sanc.violence.covar)

###plot violence effects on sanctions support using ORDINAL coding

## Difference in means plots

par(mfrow=c(1,1))
plot(effect_sanc_violence.ord,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_sanc_violence.ord_95$conf.int[1],3,diff.in.means_sanc_violence.ord_95$conf.int[2],3,lwd=2)
segments(diff.in.means_sanc_violence.ord_90$conf.int[1],3,diff.in.means_sanc_violence.ord_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Support for Sanctions",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)


##########################Appendix ####################

## EXPLANATION: LINES 1622 TO 1648 GIVES APPENDIX FIGURE 3 PLOT *lines 1567 to 1617 required first ### 

## EXPLANATION: LINES 1692 TO 1718 GIVES APPENDIX FIGURE 4 PLOT *lines 1651 to 1687 required first ### 

## EXPLANATION: LINES 1833 TO 1860 GIVES APPENDIX FIGURE 5 PLOT *lines 1720 to 1831 required first ### 

#########################Additional Analyses in Appendix for Study 2##################################

##effects of mobile treatment on blame

## Bivariate OLS (no controls)

# blame protesters outcome - mobile is not significant

reg.blame_protesters.mobile.bivar <- lm(blameprotestor ~ mobile, data=exp)
summary(reg.blame_protesters.mobile.bivar)

# blame police outcome - mobile is not significant

reg.blame_police.mobile.bivar <- lm(blamepolice ~ mobile, data=exp)
summary(reg.blame_police.mobile.bivar)

# blame government outcome - mobile is not significant

reg.blame_gov.mobile.bivar <- lm(blamegov ~ mobile, data=exp)
summary(reg.blame_gov.mobile.bivar)


## OLS with controls


# blame protesters outcome - mobile is not significant, violence is significant

reg.blame_protesters.covar <- lm(blameprotestor ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.blame_protesters.covar)

#~# blame police outcome - mobile is not significant, violence *is significant

reg.blame_police.covar <- lm(blamepolice ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.blame_police.covar)

#~# blame government outcome - mobile is not significant, violence *is significant

reg.blame_gov.covar <- lm(blamegov ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.blame_gov.covar)

## Difference in means equations

diff.in.means_mobile_blame_protesters_95 <- t.test(exp[exp$mobile==1,]$blameprotestor, exp[exp$mobile==0,]$blameprotestor)
diff.in.means_mobile_blame_protesters_90 <- t.test(exp[exp$mobile==1,]$blameprotestor, exp[exp$mobile==0,]$blameprotestor, conf.level = 0.9)
effect_mobile_blame_protesters <- diff.in.means_mobile_blame_protesters_95$estimate[1]-diff.in.means_mobile_blame_protesters_95$estimate[2]

diff.in.means_mobile_blame_police_95 <- t.test(exp[exp$mobile==1,]$blamepolice, exp[exp$mobile==0,]$blamepolice)
diff.in.means_mobile_blame_police_90 <- t.test(exp[exp$mobile==1,]$blamepolice, exp[exp$mobile==0,]$blamepolice, conf.level = 0.9)
effect_mobile_blame_police <- diff.in.means_mobile_blame_police_95$estimate[1]-diff.in.means_mobile_blame_police_95$estimate[2]

diff.in.means_mobile_blame_gov_95 <- t.test(exp[exp$mobile==1,]$blamegov, exp[exp$mobile==0,]$blamegov)
diff.in.means_mobile_blame_gov_90 <- t.test(exp[exp$mobile==1,]$blamegov, exp[exp$mobile==0,]$blamegov, conf.level = 0.9)
effect_mobile_blame_gov <- diff.in.means_mobile_blame_gov_95$estimate[1]-diff.in.means_mobile_blame_gov_95$estimate[2]


## Difference in means plots

# PLOT FOR APPENDIX FIGURE 3
#pdf("study2_mobile_blame_w90_fixed.pdf")
par(mfrow=c(3,1))
plot(effect_mobile_blame_protesters,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_protesters_95$conf.int[1],3,diff.in.means_mobile_blame_protesters_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_protesters_90$conf.int[1],3,diff.in.means_mobile_blame_protesters_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Protesters",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect_mobile_blame_police,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_police_95$conf.int[1],3,diff.in.means_mobile_blame_police_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_police_90$conf.int[1],3,diff.in.means_mobile_blame_police_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Police",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect_mobile_blame_gov,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_gov_95$conf.int[1],3,diff.in.means_mobile_blame_gov_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_gov_90$conf.int[1],3,diff.in.means_mobile_blame_gov_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Government",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
#dev.off()


##effects of violence treatment on blame


## Bivariate OLS (no controls)

# blame protesters outcome - violence is significant, p = 5.78e-06

reg.blame_protesters.violence.bivar <- lm(blameprotestor ~ violence, data=exp)
summary(reg.blame_protesters.violence.bivar)

#~# blame police outcome - violence *is significant, p = 0.0127

reg.blame_police.violence.bivar <- lm(blamepolice ~ violence, data=exp)
summary(reg.blame_police.violence.bivar)

#~# blame government outcome - violence *is significant, p = 0.0162

reg.blame_gov.violence.bivar <- lm(blamegov ~ violence, data=exp)
summary(reg.blame_gov.violence.bivar)


## OLS with controls 


## Difference in means equations

diff.in.means_violence_blame_protesters_95 <- t.test(exp[exp$violence==1,]$blameprotestor, exp[exp$violence==0,]$blameprotestor)
diff.in.means_violence_blame_protesters_90 <- t.test(exp[exp$violence==1,]$blameprotestor, exp[exp$violence==0,]$blameprotestor, conf.level = 0.9)
effect_violence_blame_protesters <- diff.in.means_violence_blame_protesters_95$estimate[1]-diff.in.means_violence_blame_protesters_95$estimate[2]

diff.in.means_violence_blame_police_95 <- t.test(exp[exp$violence==1,]$blamepolice, exp[exp$violence==0,]$blamepolice)
diff.in.means_violence_blame_police_90 <- t.test(exp[exp$violence==1,]$blamepolice, exp[exp$violence==0,]$blamepolice, conf.level = 0.9)
effect_violence_blame_police <- diff.in.means_violence_blame_police_95$estimate[1]-diff.in.means_violence_blame_police_95$estimate[2]

diff.in.means_violence_blame_gov_95 <- t.test(exp[exp$violence==1,]$blamegov, exp[exp$violence==0,]$blamegov)
diff.in.means_violence_blame_gov_90 <- t.test(exp[exp$violence==1,]$blamegov, exp[exp$violence==0,]$blamegov, conf.level = 0.9)
effect_violence_blame_gov <- diff.in.means_violence_blame_gov_95$estimate[1]-diff.in.means_violence_blame_gov_95$estimate[2]


## Difference in means plots

# PLOT FOR APPENDIX FIGURE 4
#pdf("study2_violence_blame_w90_fixed.pdf")
par(mfrow=c(3,1))
plot(effect_violence_blame_protesters,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_protesters_95$conf.int[1],3,diff.in.means_violence_blame_protesters_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_protesters_90$conf.int[1],3,diff.in.means_violence_blame_protesters_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Protesters",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect_violence_blame_police,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_police_95$conf.int[1],3,diff.in.means_violence_blame_police_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_police_90$conf.int[1],3,diff.in.means_violence_blame_police_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Police",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect_violence_blame_gov,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_gov_95$conf.int[1],3,diff.in.means_violence_blame_gov_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_gov_90$conf.int[1],3,diff.in.means_violence_blame_gov_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Government",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
#dev.off()

##effects of mobile x violence on trustworthiness of video

# Bivariate OLS, neither mobile nor violence are significant

reg.binary_trust.mobile.bivar <- lm(trustworthy.binary ~  mobile, data=exp)
summary(reg.binary_trust.mobile.bivar)

reg.binary_trust.violence.bivar <- lm(trustworthy.binary ~ violence, data=exp)
summary(reg.binary_trust.violence.bivar)

## Simple OLS with interaction term (no controls) 

# binary trustworthiness outcome - mobile, violence, interaction not significant

reg.binary_trust.violencemobile.inter <- lm(trustworthy.binary ~  mobile + violence + mobile*violence, data=exp)
summary(reg.binary_trust.violencemobile.inter)

# ordinal trust outcome - mobile, violence, interaction not significant

reg.ord_trust.violencemobile.inter <- lm(trustworthy.ordinal ~  mobile + violence + mobile*violence, data=exp)
summary(reg.ord_trust.violencemobile.inter)


## OLS with controls and interaction term

# binary trustworthiness outcome - violence is significant, p = 0.0713, but mobile and interaction not significant

reg.binary_trust.violencemobile.covar.inter <- lm(trustworthy.binary ~ mobile + violence + mobile*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.binary_trust.violencemobile.covar.inter)

# ordinal sanctions outcome - violence is significant, p = 0.06344, but mobile and interaction not significant

reg.ord_trust.violencemobile.covar.inter <- lm(trustworthy.ordinal ~ mobile + violence + mobile*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp)
summary(reg.ord_trust.violencemobile.covar.inter)

## Conditional effects plots 

# trustworthiness

# binary, marginal effect of mobile conditional on V = 1
vcov(reg.binary_trust.violencemobile.inter)
# var for mobile is 0.0006421250
# var for interaction term is 0.0012739550
# cov for mobile and interaction is -0.0006421250


margeffect_trust <-  summary(reg.binary_trust.violencemobile.inter)$coef[2] + summary(reg.binary_trust.violencemobile.inter)$coef[4]
var_margeffect_trust <- vcov(reg.binary_trust.violencemobile.inter)[2,2] + vcov(reg.binary_trust.violencemobile.inter)[4,4] + (2 * vcov(reg.binary_trust.violencemobile.inter)[4,2])
se_margeffect_trust <- sqrt(var_margeffect_trust)

ci_upper_trust <- margeffect_trust + (1.96*se_margeffect_trust)
ci_lower_trust <- margeffect_trust - (1.96*se_margeffect_trust)

ci_upper_trust90 <- margeffect_trust + (1.645*se_margeffect_trust)
ci_lower_trust90 <- margeffect_trust - (1.645*se_margeffect_trust)



# binary, marginal effect of violence conditional on M = 1
vcov(reg.binary_trust.violencemobile.inter)
# var for violence is 0.0006454418
# var for interaction term is 0.0012739550
# cov for violence and interaction is -0.0006454418 


margeffect_trust.vio <-  summary(reg.binary_trust.violencemobile.inter)$coef[3] + summary(reg.binary_trust.violencemobile.inter)$coef[4]
var_margeffect_trust.vio <- vcov(reg.binary_trust.violencemobile.inter)[3,3] + vcov(reg.binary_trust.violencemobile.inter)[4,4] + (2 * vcov(reg.binary_trust.violencemobile.inter)[4,3])
se_margeffect_trust.vio <- sqrt(var_margeffect_trust.vio)

ci_upper_trust.vio <- margeffect_trust.vio + (1.96*se_margeffect_trust.vio)
ci_lower_trust.vio <- margeffect_trust.vio - (1.96*se_margeffect_trust.vio)

ci_upper_trust90.vio <- margeffect_trust.vio + (1.645*se_margeffect_trust.vio)
ci_lower_trust90.vio <- margeffect_trust.vio - (1.645*se_margeffect_trust.vio)

# making the plots

#binary_trust.inter.mobile <- summary(reg.binary_trust.violencemobile.inter)$coefficients[2]
binary_trust.inter.violence <- summary(reg.binary_trust.violencemobile.inter)$coefficients[3]
binary_trust.inter.margeffect.vio <- margeffect_trust.vio
binary_trust.inter.margeffect <- margeffect_trust

#binary_trust.inter.mobile.err <- summary(reg.binary_trust.violencemobile.inter)$coefficients[,2][2]
#li_binary_trust.inter.mobile <- binary_trust.inter.mobile-1.96 *binary_trust.inter.mobile.err
#ui_binary_trust.inter.mobile <- binary_trust.inter.mobile+1.96 *binary_trust.inter.mobile.err

binary_trust.inter.violence.err <- summary(reg.binary_trust.violencemobile.inter)$coefficients[,2][3]
li_binary_trust.inter.violence <- binary_trust.inter.violence-1.96 *binary_trust.inter.violence.err
ui_binary_trust.inter.violence<- binary_trust.inter.violence+1.96 *binary_trust.inter.violence.err
li_binary_trust.inter.violence90 <- binary_trust.inter.violence-1.645 *binary_trust.inter.violence.err
ui_binary_trust.inter.violence90 <- binary_trust.inter.violence+1.645 *binary_trust.inter.violence.err


li_binary_trust.inter.margeffect.vio <- ci_lower_trust.vio
ui_binary_trust.inter.margeffect.vio <- ci_upper_trust.vio
li_binary_trust.inter.margeffect90.vio <- ci_lower_trust90.vio
ui_binary_trust.inter.margeffect90.vio <- ci_upper_trust90.vio


li_binary_trust.inter.margeffect <- ci_lower_trust
ui_binary_trust.inter.margeffect <- ci_upper_trust
li_binary_trust.inter.margeffect90 <- ci_lower_trust90
ui_binary_trust.inter.margeffect90 <- ci_upper_trust90

binary_trust.inter_coef <- c(binary_trust.inter.violence, binary_trust.inter.margeffect.vio, binary_trust.inter.margeffect)
binary_trust.inter_li <- c(li_binary_trust.inter.violence, li_binary_trust.inter.margeffect.vio, li_binary_trust.inter.margeffect)
binary_trust.inter_ui <- c(ui_binary_trust.inter.violence, ui_binary_trust.inter.margeffect.vio, ui_binary_trust.inter.margeffect)
binary_trust.inter_li90 <- c(li_binary_trust.inter.violence90, li_binary_trust.inter.margeffect90.vio, li_binary_trust.inter.margeffect90)
binary_trust.inter_ui90 <- c(ui_binary_trust.inter.violence90, ui_binary_trust.inter.margeffect90.vio, ui_binary_trust.inter.margeffect90)

treatment <- c("violence, m=0", "violence, m=1", "mobile, v=1")
df_binary_trust.inter <- as.data.frame(cbind(treatment, binary_trust.inter_coef, binary_trust.inter_li90, binary_trust.inter_ui90, binary_trust.inter_li, binary_trust.inter_ui))

# PLOT FOR APPENDIX FIGURE 5
#pdf("study2_trust_new_w90.pdf")
par(mfrow=c(3,1))

plot(binary_trust.inter.violence,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.violence,3,ui_binary_trust.inter.violence,3,lwd=2)
segments(li_binary_trust.inter.violence90,3,ui_binary_trust.inter.violence90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(binary_trust.inter.margeffect.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.margeffect.vio,3,ui_binary_trust.inter.margeffect.vio,3,lwd=2)
segments(li_binary_trust.inter.margeffect90.vio,3,ui_binary_trust.inter.margeffect90.vio,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(binary_trust.inter.margeffect,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.margeffect,3,ui_binary_trust.inter.margeffect,3,lwd=2)
segments(li_binary_trust.inter.margeffect90,3,ui_binary_trust.inter.margeffect90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
#dev.off()

#########################Checks##################################

## EXPLANATION: LINES 2097 TO 2132 GIVES APPENDIX FIGURE 6 PLOT *lines 2024-2070 required first ### 

## EXPLANATION: LINES 2266 TO 2317 GIVES APPENDIX FIGURE 7 PLOT *lines 2136-2234 required first ### 

## EXPLANATION: LINES 2410 TO 2463 GIVES APPENDIX FIGURE 8 PLOT *lines 2320-2378 required first ### 

## EXPLANATION: LINES 1897 TO 2007 GIVES DATA FOR APPENDIX TABLES 6 AND 7 ### 

## EXPLANATION: LINES 1881 TO 1888 GIVES APPENDIX MANIPULATION CHECK ESTIMATE OF ~90% (ON APPENDIX PAGE 7) ### 

## EXPLANATION: LINES 1889 TO 1894 GIVES APPENDIX ATTENTION CHECK ESTIMATE OF ~95% (ON APPENDIX PAGE 7) ### 

## EXPLANATION: LINES 2007 TO 2016 GIVES APPENDIX ROBUSTNESS CHECK ESTIMATE OF ~26% (ON APPENDIX PAGE 11) ### 


## MANIPULATION CHECK, ATTENTION CHECK, BALANCE CHECKS AND ROBUSTNESS CHECKS

##manipulation check
exp$mani_check <- ifelse(exp$manipulation.check == 1 & exp$mobile == 0, 1, NA)
exp$mani_check <- ifelse(exp$manipulation.check == 2 & exp$mobile == 1, 1, exp$mani_check)
exp$mani_check <- ifelse(is.na(exp$mani_check), 0, exp$mani_check)
mean(exp$mani_check) # 90% of respondents responded to the treatment

log_print("Study 2, manipulation check")
log_print(mean(exp$mani_check))
##attention check
exp$attention <- ifelse (as.numeric(as.character(exp$attention)) == 3, 1, 0)
mean(na.omit(exp$attention)) # 95% of respondents were paying attention

log_print("Study 2, attention check")
log_print(mean(na.omit(exp$attention)))
############Balance Tests#################

# balance checks on control variables

mean(exp[exp$mobile==1,]$age)
mean(exp[exp$mobile==0,]$age)
sd(exp[exp$mobile==1,]$age)
sd(exp[exp$mobile==0,]$age)
mean(exp[exp$violence==1,]$age)
mean(exp[exp$violence==0,]$age)
sd(exp[exp$violence==1,]$age)
sd(exp[exp$violence==0,]$age)

mean(exp[exp$mobile==1,]$male)
mean(exp[exp$mobile==0,]$male)
sd(exp[exp$mobile==1,]$male)
sd(exp[exp$TV==1,]$male)
mean(exp[exp$violence==1,]$male)
mean(exp[exp$violence==0,]$male)
sd(exp[exp$violence==1,]$male)
sd(exp[exp$violence==0,]$male)

mean(exp[exp$mobile==1,]$white)
mean(exp[exp$mobile==0,]$white)
sd(exp[exp$mobile==1,]$white)
sd(exp[exp$mobile==0,]$white)
mean(exp[exp$violence==1,]$white)
mean(exp[exp$violence==0,]$white)
sd(exp[exp$violence==1,]$white)
sd(exp[exp$violence==0,]$white)

mean(exp[exp$mobile==1,]$vote_presidential)
mean(exp[exp$mobile==0,]$vote_presidential)
sd(exp[exp$mobile==1,]$vote_presidential)
sd(exp[exp$mobile==0,]$vote_presidential)
mean(exp[exp$violence==1,]$vote_presidential)
mean(exp[exp$violence==0,]$vote_presidential)
sd(exp[exp$violence==1,]$vote_presidential)
sd(exp[exp$violence==0,]$vote_presidential)

mean(exp[exp$mobile==1,]$republican)
sd(exp[exp$mobile==1,]$republican)
mean(exp[exp$mobile==0,]$republican)
sd(exp[exp$mobile==0,]$republican)
mean(exp[exp$violence==1,]$republican)
sd(exp[exp$violence==1,]$republican)
mean(exp[exp$violence==0,]$republican)
sd(exp[exp$violence==0,]$republican)

mean(exp[exp$mobile==0,]$college_degree)
sd(exp[exp$mobile==0,]$college_degree)
mean(exp[exp$mobile==1,]$college_degree)
sd(exp[exp$mobile==1,]$college_degree)
mean(exp[exp$violence==0,]$college_degree)
sd(exp[exp$violence==0,]$college_degree)
mean(exp[exp$violence==1,]$college_degree)
sd(exp[exp$violence==1,]$college_degree)

mean(exp[exp$mobile==1,]$mid_income)
sd(exp[exp$mobile==1,]$mid_income)
mean(exp[exp$mobile==0,]$mid_income)
sd(exp[exp$mobile==0,]$mid_income)
mean(exp[exp$violence==1,]$mid_income)
sd(exp[exp$violence==1,]$mid_income)
mean(exp[exp$violence==0,]$mid_income)
sd(exp[exp$violence==0,]$mid_income)

mean(as.numeric(as.character(exp[exp$mobile==1,]$attention_news)))
mean(as.numeric(as.character(exp[exp$mobile==0,]$attention_news)))
sd(as.numeric(as.character(exp[exp$mobile==1,]$attention_news)))
sd(as.numeric(as.character(exp[exp$mobile==0,]$attention_news)))
mean(as.numeric(as.character(exp[exp$violence==1,]$attention_news)))
mean(as.numeric(as.character(exp[exp$violence==0,]$attention_news)))
sd(as.numeric(as.character(exp[exp$violence==1,]$attention_news)))
sd(as.numeric(as.character(exp[exp$violence==0,]$attention_news)))

mean(na.omit(exp[exp$mobile==1,]$Latin_thermometer))
mean(na.omit(exp[exp$mobile==0,]$Latin_thermometer))
sd(na.omit(exp[exp$mobile==1,]$Latin_thermometer))
sd(na.omit(exp[exp$mobile==0,]$Latin_thermometer))
mean(na.omit(exp[exp$violence==1,]$Latin_thermometer))
mean(na.omit(exp[exp$violence==0,]$Latin_thermometer))
sd(na.omit(exp[exp$violence==1,]$Latin_thermometer))
sd(na.omit(exp[exp$violence==0,]$Latin_thermometer))

mean(na.omit(exp[exp$mobile==1,]$USpolice_confidence))
sd(na.omit(exp[exp$mobile==1,]$USpolice_confidence))
mean(na.omit(exp[exp$mobile==0,]$USpolice_confidence))
sd(na.omit(exp[exp$mobile==0,]$USpolice_confidence))
mean(na.omit(exp[exp$violence==1,]$USpolice_confidence))
sd(na.omit(exp[exp$violence==1,]$USpolice_confidence))
mean(na.omit(exp[exp$violence==0,]$USpolice_confidence))
sd(na.omit(exp[exp$violence==0,]$USpolice_confidence))

mean(na.omit(exp[exp$mobile==1,]$USgovt_confidence))
sd(na.omit(exp[exp$mobile==1,]$USgovt_confidence))
mean(na.omit(exp[exp$mobile==0,]$USgovt_confidence))
sd(na.omit(exp[exp$mobile==0,]$USgovt_confidence))
mean(na.omit(exp[exp$violence==1,]$USgovt_confidence))
sd(na.omit(exp[exp$violence==1,]$USgovt_confidence))
mean(na.omit(exp[exp$violence==0,]$USgovt_confidence))
sd(na.omit(exp[exp$violence==0,]$USgovt_confidence))

mean(na.omit(exp[exp$mobile==1,]$USmedia_confidence))
sd(na.omit(exp[exp$mobile==1,]$USmedia_confidence))
mean(na.omit(exp[exp$mobile==0,]$USmedia_confidence))
sd(na.omit(exp[exp$mobile==0,]$USmedia_confidence))
mean(na.omit(exp[exp$violence==1,]$USmedia_confidence))
sd(na.omit(exp[exp$violence==1,]$USmedia_confidence))
mean(na.omit(exp[exp$violence==0,]$USmedia_confidence))
sd(na.omit(exp[exp$violence==0,]$USmedia_confidence))

# Robustness Checks - removing those who knew it was Venezuela

# finding venezuela answers to the country location check
# most venezuela answers have "venez" in them...any venezuela guesses with other spelling?
vterms <- exp$locationcheckcountry[grepl("v", exp$locationcheckcountry, ignore.case = TRUE)]
vterms[-which(grepl("venez", vterms, ignore.case = TRUE))]
vterms_minus_venez <- vterms[grepl("venez", exp$locationcheckcountry, ignore.case = FALSE)]
# veniz, venz, venuz, Venuaz, Venuez, Venaz are other spellings
exp$venezuela <- ifelse(grepl("venez", exp$locationcheckcountry, ignore.case = TRUE) | grepl("veniz", exp$locationcheckcountry, ignore.case = TRUE) | grepl("venz", exp$locationcheckcountry, ignore.case = TRUE) | grepl("venuz", exp$locationcheckcountry, ignore.case = TRUE) | grepl("venuaz", exp$locationcheckcountry, ignore.case = TRUE) | grepl("venuez", exp$locationcheckcountry, ignore.case = TRUE) | grepl("venaz", exp$locationcheckcountry, ignore.case = TRUE), 1, 0)
mean(exp$venezuela) # 26% of respondents knew this was Venezuela

log_print("Study 2 Robustness Check, aware of location of videos")
log_print(mean(exp$venezuela))

unaware <- exp[exp$venezuela == 0,]
nrow(unaware)

## Bivariate OLS (no controls)

# binary sanctions outcome - violence is still significant, p changes from 0.0036 to .0244

reg.binary_sanc.violence.bivar.rc <- lm(sanctions.binary ~ violence, data=unaware)
summary(reg.binary_sanc.violence.bivar.rc)

# ordinal sanctions outcome - violence is still significant, p changes from 0.00128 to 0.05

reg.ord_sanc.violence.bivar.rc <- lm(sanctions.ordinal ~ violence, data=unaware)
summary(reg.ord_sanc.violence.bivar.rc)

# binary learning more - violence is still not significant

reg.learn.violence.bivar.rc <- lm(learnmore ~ violence, data=unaware)
summary(reg.learn.violence.bivar.rc)



## OLS with controls

# binary sanctions outcome - violence is still significant, p changes from 0.002297 to .013861

reg.binary_sanc.violence.covar.rc <- lm(sanctions.binary ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.binary_sanc.violence.covar.rc)

# ordinal sanctions outcome - violence is still significant, p changes from 0.000601 to .03844

reg.ord_sanc.violence.covar.rc <- lm(sanctions.ordinal ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.ord_sanc.violence.covar.rc)

# binary learning more - violence is still not significant

reg.learn.violence.covar.rc <- lm(learnmore ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.learn.violence.covar.rc)



## Difference in means equations

diff.in.means_sanc_violence.rc_95 <- t.test(unaware[unaware$violence==1,]$sanctions.binary, unaware[unaware$violence==0,]$sanctions.binary)
diff.in.means_sanc_violence.rc_90 <- t.test(unaware[unaware$violence==1,]$sanctions.binary, unaware[unaware$violence==0,]$sanctions.binary, conf.level = 0.9)
effect_sanc_violence.rc <- diff.in.means_sanc_violence.rc_95$estimate[1]-diff.in.means_sanc_violence.rc_95$estimate[2]

diff.in.means_learn_violence.rc_95 <- t.test(unaware[unaware$violence==1,]$learnmore, unaware[unaware$violence==0,]$learnmore)
diff.in.means_learn_violence.rc_90 <- t.test(unaware[unaware$violence==1,]$learnmore, unaware[unaware$violence==0,]$learnmore, conf.level = 0.9)
effect_learn_violence.rc <- diff.in.means_learn_violence.rc_95$estimate[1]-diff.in.means_learn_violence.rc_95$estimate[2]



## Difference in means plots

#pdf("study2_violence_rc.pdf")
par(mfrow=c(2,1))
plot(effect_sanc_violence.rc,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_sanc_violence.rc_95$conf.int[1],3,diff.in.means_sanc_violence.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_sanc_violence.rc_90$conf.int[1],3,diff.in.means_sanc_violence.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Support for Sanctions",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect_learn_violence.rc,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_learn_violence.rc_95$conf.int[1],3,diff.in.means_learn_violence.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_learn_violence.rc_90$conf.int[1],3,diff.in.means_learn_violence.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect of Violence on Interest in Learning More",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
#dev.off()

# PLOT FOR APPENDIX FIGURE 6
#pdf("study2_violence_RCvenez_w90.pdf")
par(mfrow=c(2,2))

plot(effect_sanc_violence,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_sanc_violence_95$conf.int[1],3,diff.in.means_sanc_violence_95$conf.int[2],3,lwd=2)
segments(diff.in.means_sanc_violence_90$conf.int[1],3,diff.in.means_sanc_violence_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Sanctions Support \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.9,line=2.3)

plot(effect_sanc_violence.rc,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_sanc_violence.rc_95$conf.int[1],3,diff.in.means_sanc_violence.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_sanc_violence.rc_90$conf.int[1],3,diff.in.means_sanc_violence.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Sanctions Support \n no Venezuela recognition, N = 1293",side=3,cex=.9,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.9,line=2.3)

plot(effect_learn_violence,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_learn_violence_95$conf.int[1],3,diff.in.means_learn_violence_95$conf.int[2],3,lwd=2)
segments(diff.in.means_learn_violence_90$conf.int[1],3,diff.in.means_learn_violence_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Interest in Learning More \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.9,line=2.3)

plot(effect_learn_violence.rc,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_learn_violence.rc_95$conf.int[1],3,diff.in.means_learn_violence.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_learn_violence.rc_90$conf.int[1],3,diff.in.means_learn_violence.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("  Effect on Interest in Learning More \n no Venezuela recognition, N = 1293",side=3,cex=.9,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.9,line=2.3)
#dev.off()


### robustness check for regressions with interaction effects on sanctions, learning more


## Simple OLS with interaction term (no controls) 

# binary sanctions outcome - violence is still significant, p changes from 0.0139 to 0.0352, mobile and interaction are still not significant

reg.binary_sanc.violencemobile.inter.rc <- lm(sanctions.binary ~  mobile + violence + mobile*violence, data=unaware)
summary(reg.binary_sanc.violencemobile.inter.rc)

# ordinal sanctions outcome - violence is still significant (although at alpha = 0.1 level), p changes from 0.0114 to 0.0766, mobile and interaction still not significant

reg.ord_sanc.violencemobile.inter.rc <- lm(sanctions.ordinal ~ mobile + violence + mobile*violence, data=unaware)
summary(reg.ord_sanc.violencemobile.inter.rc)

# binary learning more - violence is still not significant, mobile is still not significant, interaction still is not significant

reg.learn.violencemobile.inter.rc <- lm(learnmore ~ mobile + violence +  mobile*violence, data=unaware)
summary(reg.learn.violencemobile.inter.rc)


## OLS with controls and interaction term

# binary sanctions outcome - violence is still significant, p changes from 0.0009502 to 0.019896

reg.binary_sanc.violencemobile.covar.inter.rc <- lm(sanctions.binary ~ mobile + violence + mobile*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.binary_sanc.violencemobile.covar.inter.rc)

# ordinal sanctions outcome - violence is still significant (albeit at alpha = 0.1 level), p changes from 0.000612 to 0.05900

reg.ord_sanc.violencemobile.covar.inter.rc <- lm(sanctions.ordinal ~ mobile + violence + mobile*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.ord_sanc.violencemobile.covar.inter.rc)

# binary learning more - violence still is not significant

reg.learn.violencemobile.covar.inter.rc <- lm(learnmore ~ mobile + violence + mobile*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.learn.violencemobile.covar.inter.rc)

# robustness check for calculating marginal effect of mobile when V = 1

margeffect_sanc <-  summary(reg.binary_sanc.violencemobile.inter.rc)$coef[2] + summary(reg.binary_sanc.violencemobile.inter.rc)$coef[4]
var_margeffect_sanc <- vcov(reg.binary_sanc.violencemobile.inter.rc)[2,2] + vcov(reg.binary_sanc.violencemobile.inter.rc)[4,4] + (2 * vcov(reg.binary_sanc.violencemobile.inter.rc)[4,2])
se_margeffect_sanc <- sqrt(var_margeffect_sanc)

ci_upper_sanc <- margeffect_sanc + (1.96*se_margeffect_sanc)
ci_lower_sanc <- margeffect_sanc - (1.96*se_margeffect_sanc)

ci_upper_sanc90 <- margeffect_sanc + (1.645*se_margeffect_sanc)
ci_lower_sanc90 <- margeffect_sanc - (1.645*se_margeffect_sanc)


# robustness check for calculating marginal effect of violence when M = 1

margeffect_sanc.vio <-  summary(reg.binary_sanc.violencemobile.inter.rc)$coef[3] + summary(reg.binary_sanc.violencemobile.inter.rc)$coef[4]
var_margeffect_sanc.vio <- vcov(reg.binary_sanc.violencemobile.inter.rc)[3,3] + vcov(reg.binary_sanc.violencemobile.inter.rc)[4,4] + (2 * vcov(reg.binary_sanc.violencemobile.inter.rc)[4,3])
se_margeffect_sanc.vio <- sqrt(var_margeffect_sanc.vio)

ci_upper_sanc.vio <- margeffect_sanc.vio + (1.96*se_margeffect_sanc.vio)
ci_lower_sanc.vio <- margeffect_sanc.vio - (1.96*se_margeffect_sanc.vio)

ci_upper_sanc90.vio <- margeffect_sanc.vio + (1.645*se_margeffect_sanc.vio)
ci_lower_sanc90.vio <- margeffect_sanc.vio - (1.645*se_margeffect_sanc.vio)

## Conditional effects plots 

# sanctions

#binary_sanc.inter.mobile.rc <- summary(reg.binary_sanc.violencemobile.inter.rc)$coefficients[2]
binary_sanc.inter.violence.rc <- summary(reg.binary_sanc.violencemobile.inter.rc)$coefficients[3]
binary_sanc.inter.margeffect.rc.vio <- margeffect_sanc.vio
binary_sanc.inter.margeffect.rc <- margeffect_sanc

#binary_sanc.inter.mobile.err.rc <- summary(reg.binary_sanc.violencemobile.inter.rc)$coefficients[,2][2]
#li_binary_sanc.inter.mobile.rc <- binary_sanc.inter.mobile.rc-1.96 *binary_sanc.inter.mobile.err.rc
#ui_binary_sanc.inter.mobile.rc <- binary_sanc.inter.mobile.rc+1.96 *binary_sanc.inter.mobile.err.rc

binary_sanc.inter.violence.err.rc <- summary(reg.binary_sanc.violencemobile.inter.rc)$coefficients[,2][3]
li_binary_sanc.inter.violence.rc <- binary_sanc.inter.violence.rc-1.96*binary_sanc.inter.violence.err.rc
ui_binary_sanc.inter.violence.rc <- binary_sanc.inter.violence.rc+1.96*binary_sanc.inter.violence.err.rc
li_binary_sanc.inter.violence.rc90 <- binary_sanc.inter.violence.rc-1.645*binary_sanc.inter.violence.err.rc
ui_binary_sanc.inter.violence.rc90 <- binary_sanc.inter.violence.rc+1.645*binary_sanc.inter.violence.err.rc


li_binary_sanc.inter.margeffect.rc.vio <- ci_lower_sanc.vio
ui_binary_sanc.inter.margeffect.rc.vio <- ci_upper_sanc.vio
li_binary_sanc.inter.margeffect.rc.vio90 <- ci_lower_sanc90.vio
ui_binary_sanc.inter.margeffect.rc.vio90 <- ci_upper_sanc90.vio

li_binary_sanc.inter.margeffect.rc <- ci_lower_sanc
ui_binary_sanc.inter.margeffect.rc <- ci_upper_sanc
li_binary_sanc.inter.margeffect.rc90 <- ci_lower_sanc90
ui_binary_sanc.inter.margeffect.rc90 <- ci_upper_sanc90


binary_sanc.inter_coef.rc <- c(binary_sanc.inter.violence.rc, binary_sanc.inter.margeffect.rc.vio, binary_sanc.inter.margeffect.rc)
binary_sanc.inter_li.rc <- c(li_binary_sanc.inter.violence.rc, li_binary_sanc.inter.margeffect.rc.vio, li_binary_sanc.inter.margeffect.rc)
binary_sanc.inter_ui.rc <- c(ui_binary_sanc.inter.violence.rc, ui_binary_sanc.inter.margeffect.rc.vio, ui_binary_sanc.inter.margeffect.rc)
binary_sanc.inter_li.rc90 <- c(li_binary_sanc.inter.violence.rc90, li_binary_sanc.inter.margeffect.rc.vio90, li_binary_sanc.inter.margeffect.rc90)
binary_sanc.inter_ui.rc90 <- c(ui_binary_sanc.inter.violence.rc90, ui_binary_sanc.inter.margeffect.rc.vio90, ui_binary_sanc.inter.margeffect.rc90)

treatment <- c("violence, m=0", "violence, m=1", "mobile, v=1")
df_binary_sanc.inter.rc <- as.data.frame(cbind(treatment, binary_sanc.inter_coef.rc, binary_sanc.inter_li.rc90, binary_sanc.inter_ui.rc90, binary_sanc.inter_li.rc, binary_sanc.inter_ui.rc))

#pdf("study2_sanctions_rc_new.pdf")
par(mfrow=c(3,1))

plot(binary_sanc.inter.violence.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_sanc.inter.violence.rc,3,ui_binary_sanc.inter.violence.rc,3,lwd=2)
segments(li_binary_sanc.inter.violence.rc90,3,ui_binary_sanc.inter.violence.rc90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(binary_sanc.inter.margeffect.rc.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_sanc.inter.margeffect.rc.vio ,3,ui_binary_sanc.inter.margeffect.rc.vio ,3,lwd=2)
segments(li_binary_sanc.inter.margeffect.rc.vio90 ,3,ui_binary_sanc.inter.margeffect.rc.vio90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(binary_sanc.inter.margeffect.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_sanc.inter.margeffect.rc ,3,ui_binary_sanc.inter.margeffect.rc ,3,lwd=2)
segments(li_binary_sanc.inter.margeffect.rc90 ,3,ui_binary_sanc.inter.margeffect.rc90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
#dev.off()

# PLOT FOR APPENDIX FIGURE 7
#pdf("study2_sanctions_RCvenez_new_w90.pdf")
par(mfrow=c(3,2))

plot(binary_sanc.inter.violence,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_sanc.inter.violence,3,ui_binary_sanc.inter.violence,3,lwd=2)
segments(li_binary_sanc.inter.violence90,3,ui_binary_sanc.inter.violence90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0 \n all data, N = 1748",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.73,line=2.3)

plot(binary_sanc.inter.violence.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_sanc.inter.violence.rc,3,ui_binary_sanc.inter.violence.rc,3,lwd=2)
segments(li_binary_sanc.inter.violence.rc90,3,ui_binary_sanc.inter.violence.rc90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0 \n no Venezuela recognition, N = 1293",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.73,line=2.3)

plot(binary_sanc.inter.margeffect.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(li_binary_sanc.inter.margeffect.vio,3,ui_binary_sanc.inter.margeffect.vio,3,lwd=2)
segments(li_binary_sanc.inter.margeffect90.vio,3,ui_binary_sanc.inter.margeffect90.vio,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1 \n all data, N = 1748",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.73,line=2.3)

plot(binary_sanc.inter.margeffect.rc.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_sanc.inter.margeffect.rc.vio ,3,ui_binary_sanc.inter.margeffect.rc.vio ,3,lwd=2)
segments(li_binary_sanc.inter.margeffect.rc.vio90 ,3,ui_binary_sanc.inter.margeffect.rc.vio90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1 \n no Venezuela recognition, N = 1293",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.73,line=2.3)

plot(binary_sanc.inter.margeffect,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_sanc.inter.margeffect ,3,ui_binary_sanc.inter.margeffect ,3,lwd=2)
segments(li_binary_sanc.inter.margeffect90 ,3,ui_binary_sanc.inter.margeffect90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1 \n all data, N = 1748",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.73,line=2.3)

plot(binary_sanc.inter.margeffect.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_sanc.inter.margeffect.rc ,3,ui_binary_sanc.inter.margeffect.rc ,3,lwd=2)
segments(li_binary_sanc.inter.margeffect.rc90 ,3,ui_binary_sanc.inter.margeffect.rc90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1 \n no Venezuela recognition, N = 1293",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.73,line=2.3)
#dev.off()

# nearly all the same, though the marginal effect of violence when M=1 loses significance in this robustness check


# marginal effect of mobile footage given V=1, robustness check, learning more

margeffect_learn <-  summary(reg.learn.violencemobile.inter.rc)$coef[2] + summary(reg.learn.violencemobile.inter.rc)$coef[4]
var_margeffect_learn <- vcov(reg.learn.violencemobile.inter.rc)[2,2] + vcov(reg.learn.violencemobile.inter.rc)[4,4] + (2 * vcov(reg.learn.violencemobile.inter.rc)[4,2])
se_margeffect_learn <- sqrt(var_margeffect_learn)

ci_upper_learn <- margeffect_learn + (1.96*se_margeffect_learn)
ci_lower_learn <- margeffect_learn - (1.96*se_margeffect_learn)

ci_upper_learn90 <- margeffect_learn + (1.645*se_margeffect_learn)
ci_lower_learn90 <- margeffect_learn - (1.645*se_margeffect_learn)

# marginal effect of violence footage given M=1, robustness check, learning more

margeffect_learn.vio <-  summary(reg.learn.violencemobile.inter.rc)$coef[3] + summary(reg.learn.violencemobile.inter.rc)$coef[4]
var_margeffect_learn.vio <- vcov(reg.learn.violencemobile.inter.rc)[3,3] + vcov(reg.learn.violencemobile.inter.rc)[4,4] + (2 * vcov(reg.learn.violencemobile.inter.rc)[4,3])
se_margeffect_learn.vio <- sqrt(var_margeffect_learn.vio)

ci_upper_learn.vio <- margeffect_learn.vio + (1.96*se_margeffect_learn.vio)
ci_lower_learn.vio <- margeffect_learn.vio - (1.96*se_margeffect_learn.vio)

ci_upper_learn90.vio <- margeffect_learn.vio + (1.645*se_margeffect_learn.vio)
ci_lower_learn90.vio <- margeffect_learn.vio - (1.645*se_margeffect_learn.vio)



#learn.inter.mobile.rc <- summary(reg.learn.violencemobile.inter.rc)$coefficients[2]
learn.inter.violence.rc <- summary(reg.learn.violencemobile.inter.rc)$coefficients[3]
learn.inter.margeffect.rc.vio <- margeffect_learn.vio
learn.inter.margeffect.rc <- margeffect_learn

#learn.inter.mobile.err.rc <- summary(reg.learn.violencemobile.inter.rc)$coefficients[,2][2]
#li_learn.inter.mobile.rc <- learn.inter.mobile.rc-1.96 *learn.inter.mobile.err.rc
#ui_learn.inter.mobile.rc <- learn.inter.mobile.rc+1.96 *learn.inter.mobile.err.rc

learn.inter.violence.err.rc <- summary(reg.learn.violencemobile.inter.rc)$coefficients[,2][3]
li_learn.inter.violence.rc <- learn.inter.violence.rc-1.96 *learn.inter.violence.err.rc
ui_learn.inter.violence.rc <- learn.inter.violence.rc+1.96 *learn.inter.violence.err.rc
li_learn.inter.violence.rc90 <- learn.inter.violence.rc-1.645 *learn.inter.violence.err.rc
ui_learn.inter.violence.rc90 <- learn.inter.violence.rc+1.645 *learn.inter.violence.err.rc

li_learn.inter.margeffect.rc.vio <- ci_lower_learn.vio
ui_learn.inter.margeffect.rc.vio <- ci_upper_learn.vio
li_learn.inter.margeffect.rc.vio90 <- ci_lower_learn90.vio
ui_learn.inter.margeffect.rc.vio90 <- ci_upper_learn90.vio

li_learn.inter.margeffect.rc <- ci_lower_learn
ui_learn.inter.margeffect.rc <- ci_upper_learn
li_learn.inter.margeffect.rc90 <- ci_lower_learn90
ui_learn.inter.margeffect.rc90 <- ci_upper_learn90

learn.inter_coef.rc <- c(learn.inter.violence.rc, learn.inter.margeffect.rc.vio, learn.inter.margeffect.rc)
learn.inter_li.rc <- c(li_learn.inter.violence.rc, li_learn.inter.margeffect.rc.vio, li_learn.inter.margeffect.rc)
learn.inter_ui.rc <- c(ui_learn.inter.violence.rc, ui_learn.inter.margeffect.rc.vio, ui_learn.inter.margeffect.rc)
learn.inter_li.rc90 <- c(li_learn.inter.violence.rc90, li_learn.inter.margeffect.rc.vio90, li_learn.inter.margeffect.rc90)
learn.inter_ui.rc90 <- c(ui_learn.inter.violence.rc90, ui_learn.inter.margeffect.rc.vio90, ui_learn.inter.margeffect.rc90)

treatment <- c("violence, m=0", "violence, m=1", "mobile, v=1")
df_learn.inter.rc <- as.data.frame(cbind(treatment, learn.inter_coef.rc, learn.inter_li.rc90, learn.inter_ui.rc90, learn.inter_li.rc, learn.inter_ui.rc))

#pdf("study2_learn_more_rc_new.pdf")
par(mfrow=c(3,1))

plot(learn.inter.violence.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_learn.inter.violence.rc,3,ui_learn.inter.violence.rc,3,lwd=2)
segments(li_learn.inter.violence.rc90,3,ui_learn.inter.violence.rc90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(learn.inter.margeffect.rc.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_learn.inter.margeffect.rc.vio ,3,ui_learn.inter.margeffect.rc.vio ,3,lwd=2)
segments(li_learn.inter.margeffect.rc.vio90 ,3,ui_learn.inter.margeffect.rc.vio90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(learn.inter.margeffect.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_learn.inter.margeffect.rc ,3,ui_learn.inter.margeffect.rc ,3,lwd=2)
segments(li_learn.inter.margeffect.rc90 ,3,ui_learn.inter.margeffect.rc90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
#dev.off()

# PLOT FOR APPENDIX FIGURE 8
#pdf("study2_learn_more_RCvenez_new_90.pdf")
par(mfrow=c(3,2))

plot(learn.inter.violence,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_learn.inter.violence,3,ui_learn.inter.violence,3,lwd=2)
segments(li_learn.inter.violence90,3,ui_learn.inter.violence90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0 \n all data, N = 1748",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.73,line=2.3)

plot(learn.inter.violence.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_learn.inter.violence.rc,3,ui_learn.inter.violence.rc,3,lwd=2)
segments(li_learn.inter.violence.rc90,3,ui_learn.inter.violence.rc90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0 \n no Venezuela recognition, N = 1293",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.73,line=2.3)


plot(learn.inter.margeffect.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(li_learn.inter.margeffect.vio,3,ui_learn.inter.margeffect.vio,3,lwd=2)
segments(li_learn.inter.margeffect90.vio,3,ui_learn.inter.margeffect90.vio,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1 \n all data, N = 1748",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.73,line=2.3)

plot(learn.inter.margeffect.rc.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_learn.inter.margeffect.rc.vio ,3,ui_learn.inter.margeffect.rc.vio ,3,lwd=2)
segments(li_learn.inter.margeffect.rc.vio90 ,3,ui_learn.inter.margeffect.rc.vio90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1 \n no Venezuela recognition, N = 1293",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.73,line=2.3)


plot(learn.inter.margeffect,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_learn.inter.margeffect ,3,ui_learn.inter.margeffect ,3,lwd=2)
segments(li_learn.inter.margeffect90 ,3,ui_learn.inter.margeffect90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1 \n all data, N = 1748",side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.73,line=2.3)

plot(learn.inter.margeffect.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_learn.inter.margeffect.rc ,3,ui_learn.inter.margeffect.rc ,3,lwd=2)
segments(li_learn.inter.margeffect.rc90 ,3,ui_learn.inter.margeffect.rc90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1 \n no Venezuela recognition, N = 1293" ,side=3,cex=0.73,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=0.73,line=2.3)
#dev.off()

# same as before


### not featured in Appendix, but let's check robustness for the "additional analyses for Study 2) 
### (robustness checks for Figures 3-5 in the appendix)


## Bivariate OLS (no controls)

# blame protesters outcome - mobile still is not significant

reg.blame_protesters.mobile.bivar.rc <- lm(blameprotestor ~ mobile, data=unaware)
summary(reg.blame_protesters.mobile.bivar.rc)

# blame police outcome - mobile still is not significant

reg.blame_police.mobile.bivar.rc <- lm(blamepolice ~ mobile, data=unaware)
summary(reg.blame_police.mobile.bivar.rc)

# blame government outcome - mobile still is not significant

reg.blame_gov.mobile.bivar.rc <- lm(blamegov ~ mobile, data=unaware)
summary(reg.blame_gov.mobile.bivar.rc)


## OLS with controls


# blame protesters outcome - mobile still is not significant, violence is still significant but now p = 0.00091

reg.blame_protesters.covar.rc <- lm(blameprotestor ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.blame_protesters.covar.rc)

#~# blame police outcome - mobile is still not significant, violence *is still significant but now p = 0.04075

reg.blame_police.covar.rc <- lm(blamepolice ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.blame_police.covar.rc)

#~# blame government outcome - mobile is still not significant, violence *newly not significant

reg.blame_gov.covar.rc <- lm(blamegov ~ mobile + violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.blame_gov.covar.rc)

## Difference in means equations

diff.in.means_mobile_blame_protesters.rc_95 <- t.test(unaware[unaware$mobile==1,]$blameprotestor, unaware[unaware$mobile==0,]$blameprotestor)
diff.in.means_mobile_blame_protesters.rc_90 <- t.test(unaware[unaware$mobile==1,]$blameprotestor, unaware[unaware$mobile==0,]$blameprotestor, conf.level = 0.9)
effect_mobile_blame_protesters.rc <- diff.in.means_mobile_blame_protesters.rc_95$estimate[1]-diff.in.means_mobile_blame_protesters.rc_95$estimate[2]

diff.in.means_mobile_blame_police.rc_95 <- t.test(unaware[unaware$mobile==1,]$blamepolice, unaware[unaware$mobile==0,]$blamepolice)
diff.in.means_mobile_blame_police.rc_90 <- t.test(unaware[unaware$mobile==1,]$blamepolice, unaware[unaware$mobile==0,]$blamepolice, conf.level = 0.9)
effect_mobile_blame_police.rc <- diff.in.means_mobile_blame_police.rc_95$estimate[1]-diff.in.means_mobile_blame_police.rc_95$estimate[2]

diff.in.means_mobile_blame_gov.rc_95 <- t.test(unaware[unaware$mobile==1,]$blamegov, unaware[unaware$mobile==0,]$blamegov)
diff.in.means_mobile_blame_gov.rc_90 <- t.test(unaware[unaware$mobile==1,]$blamegov, unaware[unaware$mobile==0,]$blamegov, conf.level = 0.9)
effect_mobile_blame_gov.rc <- diff.in.means_mobile_blame_gov.rc_95$estimate[1]-diff.in.means_mobile_blame_gov.rc_95$estimate[2]


## Difference in means plots

#pdf("study2_mobile_blame_rc.pdf")
par(mfrow=c(3,1))
plot(effect_mobile_blame_protesters.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_protesters.rc_95$conf.int[1],3,diff.in.means_mobile_blame_protesters.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_protesters.rc_90$conf.int[1],3,diff.in.means_mobile_blame_protesters.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Protesters",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect_mobile_blame_police.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_police.rc_95$conf.int[1],3,diff.in.means_mobile_blame_police.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_police.rc_90$conf.int[1],3,diff.in.means_mobile_blame_police.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Police",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect_mobile_blame_gov.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_gov.rc_95$conf.int[1],3,diff.in.means_mobile_blame_gov.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_gov.rc_90$conf.int[1],3,diff.in.means_mobile_blame_gov.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Government",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
#dev.off()

#pdf("study2_mobile_blame_RCvenez.pdf")
par(mfrow=c(3,2))

plot(effect_mobile_blame_protesters,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_protesters_95$conf.int[1],3,diff.in.means_mobile_blame_protesters_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_protesters_90$conf.int[1],3,diff.in.means_mobile_blame_protesters_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Protesters \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)

plot(effect_mobile_blame_protesters.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_protesters.rc_95$conf.int[1],3,diff.in.means_mobile_blame_protesters.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_protesters.rc_90$conf.int[1],3,diff.in.means_mobile_blame_protesters.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Protesters \n no Venezuela recognition, N = 1293",side=3,cex=0.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=0.9,line=2.3)

plot(effect_mobile_blame_police,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_police_95$conf.int[1],3,diff.in.means_mobile_blame_police_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_police_90$conf.int[1],3,diff.in.means_mobile_blame_police_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Police \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)

plot(effect_mobile_blame_police.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_police.rc_95$conf.int[1],3,diff.in.means_mobile_blame_police.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_police.rc_90$conf.int[1],3,diff.in.means_mobile_blame_police.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Police \n no Venezuela recognition, N = 1293",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)

plot(effect_mobile_blame_gov,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_gov_95$conf.int[1],3,diff.in.means_mobile_blame_gov_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_gov_90$conf.int[1],3,diff.in.means_mobile_blame_gov_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming the Government \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)

plot(effect_mobile_blame_gov.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_mobile_blame_gov.rc_95$conf.int[1],3,diff.in.means_mobile_blame_gov.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_mobile_blame_gov.rc_90$conf.int[1],3,diff.in.means_mobile_blame_gov.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming the Government \n no Venezuela recognition, N = 1293",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)
#dev.off()


### effect of violence on blame

## Bivariate OLS (no controls)

# blame protesters outcome - violence is still significant 

reg.blame_protesters.violence.bivar.rc <- lm(blameprotestor ~ violence, data=unaware)
summary(reg.blame_protesters.violence.bivar.rc)

#~# blame police outcome - violence *is still significant

reg.blame_police.violence.bivar.rc <- lm(blamepolice ~ violence, data=unaware)
summary(reg.blame_police.violence.bivar.rc)

#~# blame government outcome - violence is *newly not significant

reg.blame_gov.violence.bivar.rc <- lm(blamegov ~ violence, data=unaware)
summary(reg.blame_gov.violence.bivar.rc)


## OLS with controls 


## Difference in means equations

diff.in.means_violence_blame_protesters.rc_95 <- t.test(unaware[unaware$violence==1,]$blameprotestor, unaware[unaware$violence==0,]$blameprotestor)
diff.in.means_violence_blame_protesters.rc_90 <- t.test(unaware[unaware$violence==1,]$blameprotestor, unaware[unaware$violence==0,]$blameprotestor, conf.level = 0.9)
effect_violence_blame_protesters.rc <- diff.in.means_violence_blame_protesters.rc_95$estimate[1]-diff.in.means_violence_blame_protesters.rc_95$estimate[2]

diff.in.means_violence_blame_police.rc_95 <- t.test(unaware[unaware$violence==1,]$blamepolice, unaware[unaware$violence==0,]$blamepolice)
diff.in.means_violence_blame_police.rc_90 <- t.test(unaware[unaware$violence==1,]$blamepolice, unaware[unaware$violence==0,]$blamepolice, conf.level = 0.9)
effect_violence_blame_police.rc <- diff.in.means_violence_blame_police.rc_95$estimate[1]-diff.in.means_violence_blame_police.rc_95$estimate[2]

diff.in.means_violence_blame_gov.rc_95 <- t.test(unaware[unaware$violence==1,]$blamegov, unaware[unaware$violence==0,]$blamegov)
diff.in.means_violence_blame_gov.rc_90 <- t.test(unaware[unaware$violence==1,]$blamegov, unaware[unaware$violence==0,]$blamegov, conf.level = 0.9)
effect_violence_blame_gov.rc <- diff.in.means_violence_blame_gov.rc_95$estimate[1]-diff.in.means_violence_blame_gov.rc_95$estimate[2]


## Difference in means plots

# blaming the government loses significance with the robustness check

#pdf("study2_violence_blame_rc.pdf")
par(mfrow=c(3,1))
plot(effect_violence_blame_protesters.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_protesters.rc_95$conf.int[1],3,diff.in.means_violence_blame_protesters.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_protesters.rc_90$conf.int[1],3,diff.in.means_violence_blame_protesters.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Protesters",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect_violence_blame_police.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_police.rc_95$conf.int[1],3,diff.in.means_violence_blame_police.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_police.rc_90$conf.int[1],3,diff.in.means_violence_blame_police.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Police",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect_violence_blame_gov.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_gov.rc_95$conf.int[1],3,diff.in.means_violence_blame_gov.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_gov.rc_90$conf.int[1],3,diff.in.means_violence_blame_gov.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Government",side=3,cex=1.4,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
#dev.off()

#pdf("study2_violence_blame_RCvenez.pdf")
par(mfrow=c(3,2))

plot(effect_violence_blame_protesters,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_protesters_95$conf.int[1],3,diff.in.means_violence_blame_protesters_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_protesters_90$conf.int[1],3,diff.in.means_violence_blame_protesters_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Protesters \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)

plot(effect_violence_blame_protesters.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_protesters.rc_95$conf.int[1],3,diff.in.means_violence_blame_protesters.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_protesters.rc_90$conf.int[1],3,diff.in.means_violence_blame_protesters.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Protesters \n no Venezuela recognition, N = 1293",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)

plot(effect_violence_blame_police,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_police_95$conf.int[1],3,diff.in.means_violence_blame_police_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_police_90$conf.int[1],3,diff.in.means_violence_blame_police_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Police \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)

plot(effect_violence_blame_police.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_police.rc_95$conf.int[1],3,diff.in.means_violence_blame_police.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_police.rc_90$conf.int[1],3,diff.in.means_violence_blame_police.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming Police \n no Venezuela recognition, N = 1293",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)

plot(effect_violence_blame_gov,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_gov_95$conf.int[1],3,diff.in.means_violence_blame_gov_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_gov_90$conf.int[1],3,diff.in.means_violence_blame_gov_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming the Government \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)

plot(effect_violence_blame_gov.rc,3,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="",xaxt="n", ylim=c(0,4))
segments(diff.in.means_violence_blame_gov.rc_95$conf.int[1],3,diff.in.means_violence_blame_gov.rc_95$conf.int[2],3,lwd=2)
segments(diff.in.means_violence_blame_gov.rc_90$conf.int[1],3,diff.in.means_violence_blame_gov.rc_90$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Effect on Blaming the Government \n no Venezuela recognition, N = 1293",side=3,cex=.9,line=1)
axis(1,at=c(seq(-10,10,by=1)),labels=c("-10", "-9","-8", "-7", "-6", "-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
mtext("Change in Blame Percentage",side=1,cex=.9,line=2.3)
#dev.off()

### effects of mobile/violence on trustworthiness

# Bivariate OLS, still, neither mobile nor violence are significant

reg.binary_trust.mobile.bivar.rc <- lm(trustworthy.binary ~  mobile, data=unaware)
summary(reg.binary_trust.mobile.bivar.rc)

reg.binary_trust.violence.bivar.rc <- lm(trustworthy.binary ~ violence, data=unaware)
summary(reg.binary_trust.violence.bivar.rc)

## Simple OLS with interaction term (no controls) 

# binary trustworthiness outcome - mobile, violence, interaction still not significant

reg.binary_trust.violencemobile.inter.rc <- lm(trustworthy.binary ~  mobile + violence + mobile*violence, data=unaware)
summary(reg.binary_trust.violencemobile.inter.rc)

# ordinal trust outcome - mobile and interaction still not significant, violence BECOMES SIGNIFICANT at alpha = 0.1 level 

reg.ord_trust.violencemobile.inter.rc <- lm(trustworthy.ordinal ~  mobile + violence + mobile*violence, data=unaware)
summary(reg.ord_trust.violencemobile.inter.rc)


## OLS with controls and interaction term

# binary trustworthiness outcome - violence is not signficant THOUGH IT WAS in the original result, mobile and interaction still not significant

reg.binary_trust.violencemobile.covar.inter.rc <- lm(trustworthy.binary ~ mobile + violence + mobile*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.binary_trust.violencemobile.covar.inter.rc)

# ordinal sanctions outcome - violence is still significant, p changes from 0.06344 to 0.0766, but mobile and interaction not significant

reg.ord_trust.violencemobile.covar.inter.rc <- lm(trustworthy.ordinal ~ mobile + violence + mobile*violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=unaware)
summary(reg.ord_trust.violencemobile.covar.inter.rc)


## Conditional effects plots 

# trustworthiness, marginal effect of mobile conditional on  V = 1


margeffect_trust.rc <-  summary(reg.binary_trust.violencemobile.inter.rc)$coef[2] + summary(reg.binary_trust.violencemobile.inter.rc)$coef[4]
var_margeffect_trust.rc <- vcov(reg.binary_trust.violencemobile.inter.rc)[2,2] + vcov(reg.binary_trust.violencemobile.inter.rc)[4,4] + (2 * vcov(reg.binary_trust.violencemobile.inter.rc)[4,2])
se_margeffect_trust.rc <- sqrt(var_margeffect_trust.rc)

ci_upper_trust.rc <- margeffect_trust.rc + (1.96*se_margeffect_trust.rc)
ci_lower_trust.rc <- margeffect_trust.rc - (1.96*se_margeffect_trust.rc)

ci_upper_trust.rc90 <- margeffect_trust.rc + (1.645*se_margeffect_trust.rc)
ci_lower_trust.rc90 <- margeffect_trust.rc - (1.645*se_margeffect_trust.rc)

# trustworthiness, marginal effect of violence conditional on  M = 1


margeffect_trust.rc.vio <-  summary(reg.binary_trust.violencemobile.inter.rc)$coef[3] + summary(reg.binary_trust.violencemobile.inter.rc)$coef[4]
var_margeffect_trust.rc.vio <- vcov(reg.binary_trust.violencemobile.inter.rc)[3,3] + vcov(reg.binary_trust.violencemobile.inter.rc)[4,4] + (2 * vcov(reg.binary_trust.violencemobile.inter.rc)[4,3])
se_margeffect_trust.rc.vio <- sqrt(var_margeffect_trust.rc.vio)

ci_upper_trust.rc.vio <- margeffect_trust.rc.vio + (1.96*se_margeffect_trust.rc.vio)
ci_lower_trust.rc.vio <- margeffect_trust.rc.vio - (1.96*se_margeffect_trust.rc.vio)

ci_upper_trust.vio.rc90 <- margeffect_trust.rc.vio + (1.645*se_margeffect_trust.rc.vio)
ci_lower_trust.vio.rc90 <- margeffect_trust.rc.vio - (1.645*se_margeffect_trust.rc.vio)



# make the plots

#binary_trust.inter.mobile.rc <- summary(reg.binary_trust.violencemobile.inter.rc)$coefficients[2]
binary_trust.inter.violence.rc <- summary(reg.binary_trust.violencemobile.inter.rc)$coefficients[3]
binary_trust.inter.margeffect.rc.vio <- margeffect_trust.rc.vio
binary_trust.inter.margeffect.rc <- margeffect_trust.rc

#binary_trust.inter.mobile.err.rc <- summary(reg.binary_trust.violencemobile.inter.rc)$coefficients[,2][2]
#li_binary_trust.inter.mobile.rc <- binary_trust.inter.mobile.rc-1.96 *binary_trust.inter.mobile.err.rc
#ui_binary_trust.inter.mobile.rc <- binary_trust.inter.mobile.rc+1.96 *binary_trust.inter.mobile.err.rc


binary_trust.inter.violence.err.rc <- summary(reg.binary_trust.violencemobile.inter.rc)$coefficients[,2][3]
li_binary_trust.inter.violence.rc <- binary_trust.inter.violence.rc-1.96 *binary_trust.inter.violence.err.rc
ui_binary_trust.inter.violence.rc <- binary_trust.inter.violence.rc+1.96 *binary_trust.inter.violence.err.rc
li_binary_trust.inter.violence.rc90 <- binary_trust.inter.violence.rc-1.645 *binary_trust.inter.violence.err.rc
ui_binary_trust.inter.violence.rc90 <- binary_trust.inter.violence.rc+1.645 *binary_trust.inter.violence.err.rc

li_binary_trust.inter.margeffect.rc.vio <- ci_lower_trust.rc.vio
ui_binary_trust.inter.margeffect.rc.vio <- ci_upper_trust.rc.vio
li_binary_trust.inter.margeffect.rc.vio90 <- ci_lower_trust.vio.rc90
ui_binary_trust.inter.margeffect.rc.vio90 <- ci_upper_trust.vio.rc90


li_binary_trust.inter.margeffect.rc <- ci_lower_trust.rc
ui_binary_trust.inter.margeffect.rc <- ci_upper_trust.rc
li_binary_trust.inter.margeffect.rc90 <- ci_lower_trust.rc90
ui_binary_trust.inter.margeffect.rc90 <- ci_upper_trust.rc90






binary_trust.inter_coef.rc <- c(binary_trust.inter.violence.rc, binary_trust.inter.margeffect.rc.vio, binary_trust.inter.margeffect.rc)
binary_trust.inter_li.rc <- c(li_binary_trust.inter.violence.rc, li_binary_trust.inter.margeffect.rc.vio, li_binary_trust.inter.margeffect.rc)
binary_trust.inter_ui.rc <- c(ui_binary_trust.inter.violence.rc, ui_binary_trust.inter.margeffect.rc.vio, ui_binary_trust.inter.margeffect.rc)
binary_trust.inter_li.rc90 <- c(li_binary_trust.inter.violence.rc90, li_binary_trust.inter.margeffect.rc.vio90, li_binary_trust.inter.margeffect.rc90)
binary_trust.inter_ui.rc90 <- c(ui_binary_trust.inter.violence.rc90, ui_binary_trust.inter.margeffect.rc.vio90, ui_binary_trust.inter.margeffect.rc90)

treatment <- c("violence, m=0", "violence, m=1", "mobile, v=1")
df_binary_trust.inter.rc <- as.data.frame(cbind(treatment, binary_trust.inter_coef.rc, binary_trust.inter_li.rc90, binary_trust.inter_ui.rc90, binary_trust.inter_li.rc, binary_trust.inter_ui.rc))



#pdf("study2_trust_rc.pdf")
par(mfrow=c(3,1))

plot(binary_trust.inter.violence.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.violence.rc,3,ui_binary_trust.inter.violence.rc,3,lwd=2)
segments(li_binary_trust.inter.violence.rc90,3,ui_binary_trust.inter.violence.rc90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(binary_trust.inter.margeffect.rc.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.margeffect.rc.vio ,3,ui_binary_trust.inter.margeffect.rc.vio ,3,lwd=2)
segments(li_binary_trust.inter.margeffect.rc.vio90 ,3,ui_binary_trust.inter.margeffect.rc.vio90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(binary_trust.inter.margeffect.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.margeffect.rc ,3,ui_binary_trust.inter.margeffect.rc ,3,lwd=2)
segments(li_binary_trust.inter.margeffect.rc90 ,3,ui_binary_trust.inter.margeffect.rc90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1",side=3,cex=1.2,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
#dev.off()

#pdf("study2_trust_RCvenez.pdf")
par(mfrow=c(3,2))

plot(binary_trust.inter.violence,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.violence,3,ui_binary_trust.inter.violence,3,lwd=2)
segments(li_binary_trust.inter.violence90,3,ui_binary_trust.inter.violence90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0 \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.9,line=2.3)

plot(binary_trust.inter.violence.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.violence.rc,3,ui_binary_trust.inter.violence.rc,3,lwd=2)
segments(li_binary_trust.inter.violence.rc90,3,ui_binary_trust.inter.violence.rc90,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 0 \n no Venezuela recognition, N = 1293",side=3,cex=.9,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.9,line=2.3)

plot(binary_trust.inter.margeffect.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.margeffect.vio,3,ui_binary_trust.inter.margeffect.vio,3,lwd=2)
segments(li_binary_trust.inter.margeffect90.vio,3,ui_binary_trust.inter.margeffect90.vio,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1 \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(binary_trust.inter.margeffect.rc.vio,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.margeffect.rc.vio ,3,ui_binary_trust.inter.margeffect.rc.vio ,3,lwd=2)
segments(li_binary_trust.inter.margeffect.rc.vio90 ,3,ui_binary_trust.inter.margeffect.rc.vio90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Violent Repression \n Conditional on M = 1 \n no Venezuela recognition, N = 1293",side=3,cex=.9,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.9,line=2.3)

plot(binary_trust.inter.margeffect,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.margeffect,3,ui_binary_trust.inter.margeffect ,3,lwd=2)
segments(li_binary_trust.inter.margeffect90,3,ui_binary_trust.inter.margeffect90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1 \n all data, N = 1748",side=3,cex=.9,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.9,line=2.3)

plot(binary_trust.inter.margeffect.rc,3,cex=3,pch=20,xlim=c(-0.6,0.6),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(li_binary_trust.inter.margeffect.rc ,3,ui_binary_trust.inter.margeffect.rc ,3,lwd=2)
segments(li_binary_trust.inter.margeffect.rc90 ,3,ui_binary_trust.inter.margeffect.rc90 ,3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Mobile Footage \n Conditional on V = 1 \n no Venezuela recognition, N = 1293",side=3,cex=.9,line=1)
axis(1,at=c(seq(-0.6,0.6,by=0.1)),labels=c("-60", "-50","-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50", "60"))
mtext("Percentage Point Change",side=1,cex=.9,line=2.3)
#dev.off()

# marginal effect of mobile conditional on V = 1 on trustworthiness of video clip for robustness check (those unaware of location) becomes significant at alpha = 0.1


## CODE FOR APPENDIX TABLE 8
#Generate table for appendix
stargazer(reg.binary_sanc.violence.bivar, reg.binary_sanc.violence.covar, reg.learn.violence.bivar, reg.learn.violence.covar)

log_print("Study 2 Regression Table, code for Appendix Table 8")
log_print(stargazer(reg.binary_sanc.violence.bivar, reg.binary_sanc.violence.covar, reg.learn.violence.bivar, reg.learn.violence.covar))

## CODE FOR APPENDIX TABLE 9
#Generate table for appendix
stargazer(reg.binary_sanc.violencemobile.inter, reg.binary_sanc.violencemobile.covar.inter, reg.learn.violencemobile.inter, reg.learn.violencemobile.covar.inter)

log_print("Study 2 Regression Table, code for Appendix Table 9")
log_print(stargazer(reg.binary_sanc.violencemobile.inter, reg.binary_sanc.violencemobile.covar.inter, reg.learn.violencemobile.inter, reg.learn.violencemobile.covar.inter))

## EXPLANATION: LINES 2934 TO 2947 GIVE APPENDIX RESULTS DISCUSSED IN SECTION K.2

# numbers reported on Appendix page 10, section K.2
reg.binary_trust.mobileonly.covar <- lm(trustworthy.binary ~ violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp[exp$mobile ==1,])
summary(reg.binary_trust.mobileonly.covar)

reg.binary_trust.tvonly.covar <- lm(trustworthy.binary ~ violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp[exp$mobile ==0,])
summary(reg.binary_trust.tvonly.covar)

# ordinal results, discussed in Appendix footnote 5

reg.ordinal_trust.mobileonly.covar <- lm(trustworthy.ordinal ~ violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp[exp$mobile ==1,])
summary(reg.ordinal_trust.mobileonly.covar)

reg.ordinal_trust.tvonly.covar <- lm(trustworthy.ordinal ~ violence + republican + vote_presidential + college_degree + as.numeric(attention_news) + male + age + white + mid_income + Latin_thermometer + USpolice_confidence + USgovt_confidence + USmedia_confidence, data=exp[exp$mobile ==0,])
summary(reg.ordinal_trust.tvonly.covar)

log_close()

sink('TandSC-log.txt', append=TRUE)

# View results
writeLines(readLines(logfile))

sink()
